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Abstract 


The  complex  response  of  a  granular  material  to  loading  and  unloading  results 
directly  from  its  particulate  nature.  The  independent  particle  motions  affect 
load  transferal  among  neighboring  grains  and  alters  the  microstructure  of  the 
material.  Constituent  particles  typically  transfer  load  through  shear  and  normal 
forces  at  contacts  with  neighboring  particles,  causing  the  overall  stress  to  be 
unevenly  distributed  in  the  material.  This  research  focuses  on  the  development  of 
a  mathematical  model  to  derive  the  response  of  a  macroscopic  point  in  a  granular 
material  from  the  overall  response  of  collection  of  particles. 

The  mathematical  formulation  describes  a  methodology  for  deriving  the 
macroscopic  constitutive  law  for  a  granular  material  directly  from  the  discrete 
particle-to-particle  interactions.  The  important  feature  of  the  model  lies  in  the 
implicit  generation  of  the  micromechanical  response  of  a  particle  assembly  quasi- 
statically  and  not  dynamically.  The  quasi-static  approach  to  generating  the 
micromechanical  response  is  more  appropriate  for  a  broader  class  of  engineering 
applications  where  inertia  effects  are  negligible  such  as  in  creep  deformation. 
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Chapter  1 

Introduction 


1.1.  Problem  Statement 

Many  problems  of  special  interest  to  engineers  and  earth  scientists  involve 
characterization  of  the  engineering  behavior  of  granular  materials  such  as  soils. 
Granular  materials  respond  quite  differently  from  polycrystalline  materials  such  as 
metals  due  to  their  particulate  nature  and  more  random  microstructure.  Since  the 
particle-to-particle  interaction  is  frictional  in  nature,  contacts  in  granular  media  can 
easily  form  or  break  as  the  particles  move  toward  their  new  equilibrium  positions. 
From  a  modeling  standpoint,  the  prediction  of  the  new  microstructure  of  the 
material  represents  a  major  source  of  numerical  difficulty.  Consequently,  the  effects 
of  microstructure  on  the  macromechanical  response  of  granular  media  remain  poorly 
understood. 

The  distinguishing  features  of  a  granular  material’s  response  to  loading  and 
unloading  directly  result  from  the  material’s  particulate  nature.  Constituent 
particles  typically  interact  with  distinct  neighboring  granules  at  points  of  contact, 
and  during  loading/unloading,  unbalanced  shear  and  normal  contact  forces  impel 
particles  toward  their  new  equilibrium  positions.  As  a  result,  the  overall  stress 
field  becomes  unevenly  distributed  throughout  the  material.  Many  experimental, 
analytical  and  numerical  models  have  been  developed  to  study  such  particle  behavior 
and  its  effects  on  the  overall  response  of  a  granular  material. 


Several  common  experimental  methods  include  microscopic  evaluation  of  soil 
samples  and  the  loading/unloading  of  metal,  photoelastic  disk  assemblies.  Oda  [1-3] 
used  a  hardening  resin  to  fix  a  sand  sample’s  microstructure  during  deformation. 
They  reported  spatial  distributions  of  the  average  number  of  contacts  per  particle, 
orientation  of  the  long  axes  of  the  grains,  and  the  orientation  of  the  contact 
planes  between  particles.  Drescher  and  De  Josselin  de  Jong  [4]  and  Oda  et  al.  [5- 
6]  photographed  assemblies  of  photoelastic  rods  at  various  stages  of  loading  to 
determine  contact  forces,  displacements,  and  rotations  of  individual  disks.  However, 
such  analyses  are  time  consuming  and  lack  flexibility  to  run  multiple  tests  on 
precisely  identical  specimens  for  studying  the  effect  of  any  parameters. 

Analytical  models  provide  expressions  to  predict  the  nonlinear  behavior  of 
uniformly  packed  and  sized  spheres.  Deresiewiez  [7]  proposed  a  model  for  cubic 
arrays  of  uniform  spheres  that  accommodated  for  the  ultimate  failure  of  the 
assembly  and  predicted  a  hysteretic  stress-strain  behavior.  Thorton  [8]  developed  a 
general  solution  for  the  strength  of  face- centered  cubic  array  of  uniform  spheres  and 
identified  a  variety  of  failure  envelopes.  This  regular  packing  of  uniform  spheres  and 
a  few  simple  loading/unloading  paths  restrict  the  analytic  approach  to  calculating 
stresses  and  displacements  in  a  granular  media. 

In  contrast  to  experimental  and  analytical  models,  numerical  models  allow  the 
observation  of  minute  changes  in  an  idealized  material’s  structure  for  any  loading  and 
unloading  path.  Numerical  simulation  provides  access  to  detailed  micromechanical 
statistics,  particle  motions,  and  interparticle  forces.  It  also  allows  flexibility  to  alter 
loading  paths,  particle  size  distributions,  and  the  physical  properties  of  the  particles. 

Serrano  and  Rodriguez-Ortiz  [9],  Rodriguez-Ortiz  [10],  and  Kishino  [19] 
developed  a  numerical  model  for  assemblies  of  circular  disks  by  solving  for  an 
assembly  of  particle’s  deformation  based  on  governing  equations  for  each  particle  as 
it  interacts  with  its  surrounding  particles.  These  models  allow  the  specification  of 
either  a  quasi-static  stress  or  strain  history  for  the  assembly  boundary.  Serrano  and 
Rodriguez-Ortiz  calculate  particle  equilibrium  forces  and  displacements  by  assuming 
that  the  incremental  displacements  of  each  particle  center  determines  increments  of 
contact  force.  They  continually  solve  the  linear  system  of  equations  relating  the 


particle  contact  forces  to  the  particle  displacements  and  rotations.  The  contact 
theories  of  Mindlin  [11]  and  Mindlin  and  Deresiewiez  [12]  form  the  basis  of  the 
matrix  of  contact  compliances  which  must  be  updated  for  each  iteration  according 
to  the  loading  path  followed  for  each  contact.  Kishino’s  model,  on  the  other  hand, 
iteratively  displaces  and  rotates  each  particle  according  to  the  contact  stiffness 
matrix  defined  by  the  locations  of  neighboring  particles.  The  contact  stiffness  matrix 
consists  of  simply  rotating  a  constant  diagonal  matrix  of  normal  and  tangential 
contact  compliances  by  a  transformation  matrix.  The  orientation  of  the  vector 
connecting  the  centroids  of  contacting  particles  defines  this  matrix  for  each  iteration. 
The  model  applies  a  strain  history  via  the  movement  of  four  straight  walls  by 
prescribed  incremental  displacements  and  applies  a  stress  history  by  moving  the 
walls  to  attain  the  prescribed  external  stress  according  to  the  contact  stiffnesses  of 
the  boundary. 

Cundall  and  Strack  [13]  proposed  a  computationally  less  intensive  algorithm 
for  assemblies  of  particles  called  the  distinct  element  method  (DEM).  The  method 
considers  the  interaction  of  particles  as  a  transient  problem  where  states  of 
equilibrium  develop  whenever  internal  forces  balance.  The  crucial  feature  of  the 
DEM  lies  in  the  choice  of  a  small  enough  time  step  over  which  particle  velocities  and 
accelerations  may  be  assumed  constant  and  particle  disturbances  do  not  propagate 
further  than  a  disk’s  contiguous  neighbors.  Many  improvements  on  the  original 
DEM  algorithm  have  since  been  made  [14-19].  The  DEM,  however,  contains  several 
shortcomings.  The  DEM  simply  employs  a  particle’s  acceleration  to  predict  its 
direction  of  movement  when  in  fact  neighboring  particles  constrain  this  movement  to 
be  not  co-axial  with  the  particle’s  acceleration  vector.  Also,  since  the  fundamental 
idea  of  the  model  rests  on  explicit  integration  of  Newton’s  second  law  for  each 
particle,  the  model  treats  all  problems,  even  quasi-static  ones,  dynamically.  Even 
with  an  extremely  slow  rate  of  loading,  the  DEM  uses  inertial  effects  to  predict  the 
future  position  of  each  particle.  The  particle  displacements  are  calculated  based 
on  the  “bounces”  that  particles  make  against  contacting  particles,  and  not  on  the 
rotational  and  “scraping”  effects  which  dominate  a  granular  material’s  response  to 
shearing  deformations.  For  a  problem  with  a  very  slow  loading  rate  or  even  a  static 
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problem,  it  becomes  important  to  regulate  the  time  effect  in  order  for  the  inertial 
effects  to  be  negligible.  Furthermore,  the  explicit  nature  of  the  numerical  integration 
algorithm  makes  the  DEM  conditionally  stable  and  restricts  the  step  size  by  the 
natural  period  of  the  smallest  particle.  Thus,  for  even  the  simplest  quasi-static 
stress  histories,  inertial  effects  and  numerical  stability  require  many  computational 
cycles  (see,  e.g.,  [20-22]). 

In  a  different  approach  to  numerical  modeling  of  granular  materials, 
micromechanical  models  describe  the  deformation  behavior  of  an  assembly  of 
particles  using  the  continuum  concepts  of  stress  and  strain.  This  type  of 
model  derives  the  constitutive  relationship  for  a  granular  assembly  based  on  the 
microstructure  and  interaction  between  particles.  Nemat-Nassar  and  Mehrabadi  [23] 
describe  the  overall  mechanical  response  of  the  granular  mass  with  a  simple 
micromechanical  model  that  treats  the  translation  and  rotation  of  discrete  particles 
as  continuous  fields.  They  transform  a  discrete  system  into  an  equivalent  continuum 
where  continuum  concepts  of  stress  and  strain  describe  the  material’s  overall 
behavior.  The  model  decomposes  the  local  (micro)  deformation  rate  and  spin 
into  a  ‘plastic’  or  ‘inelastic’  part  which  leaves  the  soil  fabric  unchanged,  and  an 
accommodating  part  which  changes  the  soil  fabric  and  results  in  a  change  in  the 
overall  stress.  At  the  macroscopic  level,  the  notion  of  self-consistency  as  described 
by  Hill  [24,25]  relates  the  overall  nominal  stress  rate  to  the  overall  velocity  gradient 
in  terms  of  local  variables. 

Similarly,  Chang  et  al.  [26-28]  adopt  a  self-consistent  method  to  account  for 
the  effect  of  particle  separation,  particle  sliding  associated  with  large  deformations, 
and  a  nonuniform  strain  field.  Their  model  proposes  a  decomposition  of  the  general 
nonlinear  particle  displacement  field  into  a  linear  displacement  field  and  a  field 
describing  the  particle’s  movement  to  achieve  equilibrium  after  the  imposition  of  the 
linear  field.  In  addition,  the  magnitude  of  the  local  (micro)  strain  depends  largely 
on  the  local  stiffness  rather  than  on  the  location  of  a  particle  and  its  contacting 
neighbors.  In  effect,  the  macro-behavior  of  the  granular  assembly  can  be  traced  back 
to  the  micromechanical  behavior  of  a  contact.  This  model,  however,  only  subjects 
the  boundaries  of  the  representative  unit  consisting  of  many  randomly  arranged 


particles  to  displacements  compatible  with  a  uniform  strain  and  does  not  apply  a 
stress  field. 

The  modeling  approach  pursued  in  this  study  differs  somewhat  from  those 
proposed  in  previous  studies  in  that  the  nonuniform  and  locally  discontinuous 
particle  motions  directly  derive  the  macroscopic  constitutive  equations.  In  effect,  the 
model  treats  the  assembly  of  particles  as  a  separate  structure  whose  overall  response 
reflects  the  character  of  a  macroscopic  point  experiencing  uniform  deformation.  Note 
that  the  macro-element  is  interpreted  herein  as  a  homogeneously  deforming  element 
and  not  as  a  soil  sample  experiencing  non-uniform  deformation.  Viewing  the  macro¬ 
element  in  this  fashion  renders  the  resulting  model  perfectly  amenable  to  finite 
element  implementation  since  the  overall  response  of  the  particle  assembly  may  now 
be  interpreted  as  the  sampled  macroscopic  response  at  the  Gauss  integration  points. 
It  also  allows  the  representation  of  the  nonlinear  elasto-plastic  behavior  of  granular 
materials  using  laws  derived  from  particle-to-particle  contacts. 

The  central  features  of  the  model  can  be  described  as  follows.  The  model’s 
foundation  rests  on  the  conventional  hypothesis  in  computational  plasticity  by 
assuming  a  deformation-driven  problem,  i.e.,  a  given  overall  displacement  gradient 
determines  the  overall  stress  response.  The  assembly  of  particles  representing  the 
macroscopic  point  experiences  a  prescribed  overall  uniform  motion,  after  which 
the  particles  are  allowed  to  move  in  the  microscopic  sense.  During  the  particle 
motion,  contacts  are  allowed  to  form  or  break,  thus  altering  the  microstructure  of 
the  material.  The  overall  stresses  are  then  evaluated  from  the  contact  forces  that 
develop  between  the  particles  using  the  relationships  derived  from  a  virtual  work 
formulation  presented  in  [23,24,29]. 

When  the  particle  assembly  experiences  a  macroscopic  stress  history,  the 
model  iteratively  seeks  the  assembly’s  overall  uniform  motion  which  creates  contact 
forces  exactly  balancing  the  applied  stress.  The  basis  of  the  model  lies  on  the 
analytical  description  of  the  change  in  each  particle’s  displacement  (a  microscopic 
quantity)  with  respect  to  the  overall  motion  of  the  particle  £issembly  (a  macroscopic 
quantity).  In  the  terminology  of  self-consistent  models,  this  matrix  describes 
the  so-called  ‘concentration  tensor.’  To  formalize  this  complex  microscopic- 


macroscopic  relationship,  a  decomposition  similar  to  Chang  [26-28]  is  employed 
for  the  constituent  particle  motions.  Each  particle’s  movement  consists  of  a  uniform 
displacement  resulting  from  the  overall  strain  field  and  a  nonuniform,  locally 
discontinuous  motion  necessary  to  balance  the  contact  forces  among  neighboring 
particles  and  the  applied  loads.  The  exact  description  of  the  concentration  tensor, 
however,  does  not  lend  itself  to  elegant  or  efiScient  computer  implementation.  To 
circumvent  this  difficulty,  this  tensor  is  defined  with  a  secant  approximation  that 
retains  the  essential  properties  of  the  exact  description.  The  formulation  allows 
finite  motions  of  the  individual  particles,  but  the  macroscopic  material  response  is 
formulated  based  on  infinitesimal  theory.  The  model  describes  only  the  material 
overall  response  and  leaves  the  finite  rotation  and/or  large  deformation  of  the 
assembly  of  particles  for  a  future  application. 

Throughout  the  application  of  an  overall  stress  history,  the  onset  of  localization 
in  the  particle  assembly  is  investigated.  The  model’s  macroscopic  constitutive 
relation  falls  neatly  within  previous  localization  analyses  of  Rice  [30]  and  Rudnicki 
and  Rice  [31].  Localization  is  interpreted  as  an  instability  in  the  macroscopic 
constitutive  description  of  the  material’s  deformation.  Thus,  it  is  possible  for 
a  material  to  undergo  localized  deformation  even  before  the  constitutive  relation 
ceases  to  become  invertible,  or  even  before  some  related  indication  of  ideally 
plastic  response  is  met.  The  macroscopic  model’s  framework  allows  the  natural 
incorporation  of  this  localization  analysis. 

The  next  section  describes  the  format  of  the  presentation  of  the 
micromechanical-macromechanical  model. 


1.2.  Structure  of  Presentation 

Chapter  2  presents  the  kinematics  of  particle  motion  from  a  finite-displacement 
viewpoint  that  includes  particle  rotation.  This  chapter  also  details  the  contact 
model  constitutive  relation  for  describing  the  particle-to-particle  interaction.  These 
equations  will  then  be  specialized  to  a  two-dimensional  assembly  consisting  of  rigid 
circular  disks. 


Chapter  3  describes  the  field  equations  essential  for  deriving  the 
micromechanical-macromechanical  connections  between  the  overall  response  of  the 
granular  material  and  the  constituent  particle  behavior.  This  chapter  also  presents 
the  macroscopic  constitutive  relation  used  for  modeling  the  overall  response  of  the 
material  and  for  predicting  the  onset  of  localization  in  the  control  volume. 

Chapter  4  and  Chapter  5  present  the  numerical  algorithms  for  determining 
the  displacement  vector  of  each  particle  when  the  imposed  boundary  conditions  are 
strain-controlled  and  stress-controlled,  respectively.  The  notion  of  a  repeating  cell 
essential  to  the  numerical  implementation  of  the  model  will  also  be  described  as 
part  of  Chapter  4.  Chapter  5  pays  special  attention  to  the  stress  driven  problem 
and  the  exact  steps  necessary  to  solve  this  problem. 

Chapter  4  and  Chapter  5  end  with  two-dimensional  plane-strain  numerical 
simulations  on  regular  and  random  initial  packing  of  circular  disks.  The  examples 
demonstrate  the  model’s  ability  to  capture  the  periodicity  of  the  control  volume, 
its  invariance  under  rigid  body  motion,  and  its  convergence  characteristics  for  both 
boundary  conditions.  Experiments  which  include  a  localization  analysis  will  also  be 
presented. 

Chapter  6  focuses  on  model  validation.  This  chapter  includes  a  secant 
approximation  to  the  Hertzian  contact  theory  to  describe  the  contact  law  between 
two  particles.  The  performance  of  the  numerical  model  is  then  compared  with 
experimental  plane  strain  tests  on  dry  sands. 

Chapter  7  extends  the  dry  model  to  the  case  of  saturated  granular  assemblies. 
A  continuum  methodology  is  pursued  to  predict  the  pore  pressure  changes  in  fully 
saturated  granular  assemblies  as  a  result  of  changes  in  the  soil’s  microstructure. 

Chapter  8  summarizes  the  work  done.  It  also  contains  recommendations  for 
further  study  and  for  extension  of  the  model  developed  in  this  thesis. 


Chapter  2 

Micromechanical  Relations 


2.1.  Introduction 

The  fundamental  elements  of  the  model  framework  rest  on  the  interpretation 
that  the  overall  response  of  an  assembly  of  particles  represents  the  behavior  of  a 
uniformly  deforming  macroscopic  point  in  a  granular  material.  During  loading  and 
unloading  of  the  assembly,  constituent  particles  experience  locally  discontinuous 
motions  as  unbalanced  particle  contact  forces  drive  them  towards  new  equilibrium 
positions.  Within  the  assembly,  each  discrete  particle  can  be  visualized  as  a  ‘node’ 
connected  to  contiguous  nodes  via  ‘elements’  having  normal  and  tangential  degrees 
of  freedom.  When  a  contact  forms  between  two  particles,  an  element  is  created 
and  its  loading  history  tracked  until  the  particles  separate  causing  the  element 
to  disappear.  Thus,  throughout  the  loading  history,  the  number  of  nodes  in 
the  assembly  control  volume  remains  constant  but  the  number  and  orientation  of 
elements  may  fluctuate. 

Attention  in  this  chapter  focuses  on  only  the  interaction  of  two  contacting 
particles  (nodes)  within  the  assembly  control  volume  and  their  contact  element. 
First,  the  kinematics  of  the  particles  motion  are  considered  assuming  a  finite 
displacement  of  the  two  particles.  Next,  the  constitutive  equation  describing 
the  contact  element’s  elasto-plastic  behavior  will  be  presented  in  the  context  of 
plasticity  theory.  Finally,  the  constitutive  relation  will  be  integrated  assuming  a 
finite  incremental  motion. 


2.2.  Particle  kinematics 


Consider  a  finite-sized,  spherical  particle  A.  Within  particle  A,  let 
represent  a  material  vector  in  the  undeformed  configuration,  and  dx^  denote  the 
position  occupied  by  dX^  in  the  deformed  configuration.  Assume  dX'^  and  dx'^ 
are  measured  with  respect  to  the  particle  centroids  and  x^  respectively.  Thus, 
the  position  vectors  in  the  undeformed  and  deformed  configurations  of  any  point  in 
particle  A  are  given  by 

X^  =  X^  +  dX^,i  (2.1) 

x^  =  x^  +  dx^ .  (2.2) 

With  reference  to  the  undeformed  configuration,  the  deformation  gradient  F 
relates  the  material  vector  dX^  to  the  vector  dx'^  via 

dx^  =  •  dX^  ,  (2.3) 

where 

F^  =  U=V  R^,  (2.4) 

and  U  and  V  are  the  right  and  left  stretch  tensors,  respectively,  and  R^  is  the 
orthogonal  rotation  tensor.  For  a  rigid  particle,  U  =  V  =  I,  the  identity  tensor. 
The  material  vector  dX^  only  experience  a  rigid  body  rotation  and  displacement 
but  no  stretch  as  seen  by  reducing  (2.3)  to 

dx^  =  R^  •  dX^  .  (2.5) 

Thus,  particle  A  only  undergoes  rigid  body  rotation  and  translation. 

To  quantify  this  result,  let  the  centroidal  particle  displacement  d/^  be  defined 
by  dA  =  x^  —  X^.  In  the  deformed  configuration,  the  position  vector  of  any  point 
in  particle  A  then  becomes 

x^  =  X^  A  d^  A  R^-  dX^  .  (2.6) 

The  position  vector  x'^  can  be  found  by  successive  application  of  the  translation 
and  a  rigid  body  rotation  R^  of  the  material  vector  dX"^. 


9 


Now,  consider  a  contact  point  between  two  spherical  particles  A  and  B.  With 
respect  to  the  undeformed  configuration,  the  contact  point  at  C  has  a  position  vector 
given  by 

=  X^  +  =  X^  +  ,  (2.7) 

where  and  denote  the  particle  radii,  and  and  represent  the  unit 
outward  vectors  through  C.  Since  =  — iV^,  these  vectors  can  be  written  as 

j^A  ^  ^  _xAy(^^Aj^  ^B^  _  ^2.8) 

With  respect  to  the  deformed  configuration,  the  contact  point  at  C  has  a  position 
vector  given  by  (2.6)  as 

=  x^  +  .  (2.9) 

Likewise,  a  new  contact  point  c  in  the  deformed  configuration  has  a  position 
vector  given  by 

-{■  ,  (2.10) 

where  the  unit  outward  vectors  “nA  and  are  found  to  be 

nA  =  —n^  =  (x^  —  iA)l{r^  +  r^) .  (2-11) 

In  general,  the  new  contact  point  c  does  not  coincide  with  the  old  contact  point  C. 
Taking  the  magnitudes  of  equations  (2.7)  and  (2.11)  yields 

r^  +  r^=  l|i:^-±^||  =  ||a;^-a:^||.  (2.12) 

If  (2.12)  is  not  satisfied,  the  particles  will  either  separate  or  overlap.  The  initial  and 
final  indentation  becomes,  respectively, 

A  =  r^  +  rf -\\X^ -X^W  (2.13) 

6  =  r^  +  r^ -\\x^ -x^  (2.14) 


where  A,6  >  0. 

Now,  consider  the  relative  rotation  of  the  contact  point  on  particle  A  from  C 


to  c.  This  rotation  is  given  by  the  vector 


gAC  ^  qAC  gAC  ^  ^2.15) 

where 

=  {B^  ■  N^)  X  n^/  sin  \e^^\ ,  (2.16) 

and  6“^^  represents  the  magnitude  of  the  relative  rotation,  and  6^^  denotes  the 
unit  normal  to  the  plane  on  which  the  rotation  occurs.  Similar  results  are  obtained 
if  the  contact  point  is  interpret^  as  belonging  to  particle  B. 

If  particle  B  is  fixed  against  rotation  and  displacement,  the  magnitude  of  the 
contact  slip  is  given  by 

7  =  (2.17) 

which  represents  the  length  of  an  arc  on  the  surface  of  a  sphere  of  radius 
subtended  by  the  angle  0^^,  and  lying  on  a  plane  with  unit  normal  .  If  particle 
B  also  rotates,  the  slip  becomes 

'y  =  0^^ 'y  =  ll^ll .  (2.18) 


If  0^^r^  =  then  7  =  0  and  particles  A  and  B  simply  undergo  pure  rolling. 

The  rotation  vector  relationships  collapse  neatly  for  2D-plane  strain  analysis. 
Consider  the  same  particles  A  and  B,  now  represented  by  circular  disks  with  radii 
and  on  the  a:i ,  X2-plane,  respectively.  In  this  case,  the  I’otation  of  particle  A 
defined  in  (2.5)  becomes 


cos  9^  —  sin  9^ 

sin  9‘^  cos  9-^  ’ 


(2.19) 


where  9-^  denotes  the  rigid  body  rotation  of  particle  A. 

The  relative  motion  of  the  contact  points  is  now  given  by  the  scalar  angular 
rotations  9^^  and  9^^  which  from  (2.15)  and  (2.16)  can  be  written  as 


sin  9^^  =  63  •  {R^  •  N^)  X  , 


(2.20) 


Figure  1.  Slip  between  two  contacting  circular  disks. 


where  63  is  the  unit  basis  vector  in  the  X3  direction.  Expanding  the  right-hand  side 
of  (2.20)  gives 

sin  0"^^  =  cos  63  •  JV-^  y.  nA  —  sin  9^ 

=  cos  9^  sin  9^  —  sin  9^  cos  9^ 

=  sin(0*^  -  9^) ,  (2.21) 

where  9^  denotes  the  rotation  of  the  unit  normal  vector  (assumed  positive 
counterclockwise).  A  similar  result  can  be  obtained  for  particle  B. 

Thus,  equation  (2.21)  implies  that  the  relative  rotation  of  the  contact  points 
equals  the  change  in  direction  of  the  unit  normal  minus  the  total  rotation  of  the 
particle 


9AC  _qA 

0BC^qC  _qB 


(2.22) 

(2.23) 


The  particle  slip  can  be  defined  analogously  to  (2.18)  as 


j  =  9 


AC  A 


+  9 


bcb 

r 


(2.24) 


assuming  that  A,  ^  >  0.  Physically,  the  slip  7  represents  the  tangential  stretching 
of  the  contact  spring  which,  for  large  rotations,  wraps  around  the  particle  sides  as 
shown  in  Figure  1. 
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2.3.  Constitutive  relation  for  particle  contacts 


In  the  following,  Mohr- Coulomb  frictional  resistance  to  slip  governs  the 
microstructural  rearrangement  mechanism  between  two  particles.  Slip  between  two 
particles  cannot  occur  until  the  contact  shear  force  reaches  a  critical  value.  This 
particle  contact  model  is  analogous  to  crystallographic  slip  where  slip  only  occurs 
on  well  defined  slip  systems  when  a  critical  resolved  shear  stress  is  achieved.  The 
slip  plane  and  slip  direction  uniquely  define  the  slip  system. 


In  crystal  slip,  an  a  priori  known  set  of  slip  systems  constrain  the  direction  and 
location  of  slip.  The  central  task  lies  in  choosing  a  set  of  independent  active  slip 
systems  from  a  pool  of  known  linearly  dependent  potentially  active  slip  systems. 
However,  particulate  slip  occurs  in  a  direction  randomly  constrained  by  neighboring 
particles  but  on  a  known  contact  plane  with  a  normal  vector  parallel  to  the 
branch  vector  connecting  contiguous  particle  centroids.  Therefore,  in  contrast  to 
crystallographic  slip,  as  many  potential  slip  systems  exist  as  there  are  admissible 
slip  directions  —  an  infinite  set.  Since  these  systems  all  lie  in  the  same  plane  with 
only  different  slip  directions,  the  problem  of  dependent  slip  systems  does  not  exist. 
For  more  information  on  the  exact  selection  of  the  active  crystal  slip  systems  one 
may  see  Borja  and  Wren  [32]. 


In  the  following,  only  two  dimensional  circular  disk  particles  are  considered. 
The  formalism  of  plasticity  theory  has  been  adopted  for  describing  the  constitutive 
model  for  an  ideally  plastic  and  work-hardening  particle  contact.  At  the  contact 
points,  the  particles  will  exert  a  normal  force  and  a  tangential  force  /t  through 
their  contact  element.  Each  contact  element  is  represented  as  a  pair  of  springs 
with  stiffnesses  and  kT-  The  associated  deformations  are  given  by  a  normal 
indentation  6  of  the  normal  spring  and  a  slip  7  in  the  tangential  direction.  Figure  2 
shows  the  sign  conventions  for  the  forces  and  the  deformations.  In  condensed  form, 
the  contact  element  forces,  deformations,  and  stiffness  are 

fcjy  0 


JC^  = 


0  ki 


(2.25) 


If  the  spring  stiffnesses  are  constant  and  kff  ^  0,  then  the  indentation  becomes  very 
small,  and  the  particles  may  be  considered  rigid.  A  more  elaborate  constitutive 
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Figure  2.  Constitutive’' models  for  indentation  and  slip  at 
contact  points. 


model  for  the  contact  forces  is  provided  by  a  class  of  nonlinear  elastic  contact 
theories  where  the  spring  stiffnesses  are  not  constant  (see  e.g.,  [11,  12]),  but  this 
could  engender  additional  numerical  complexity. 

On  a  particle  level  the  rate-constitutive  equation  for  the  contact  element  is 
given  by 

T  =  E{S)1C  ■  if  =  H{S)V  ■  {V  -  VP) ,  (2.26) 

where  2?®  and  VP  are  the  elastic  and  plastic  components  of  V,  respectively,  and 
H{S)  is  the  Heaviside  function  which  takes  on  the  value  of  unity  when  ^  >  0  and 
zero  otherwise.  Note  that  setting  large  spring  stiffnesses  results  in  a  rigid-plastic 
contact  model. 

Plastic  tangential  slip  takes  place  when  the  tangential  contact  force  reaches 
a  critical  maximum  value  of  the  Mohr- Coulomb  frictional  resistance  given  by  the 


expression 


F=|/T|-a  +  /ivtan(^  =  0,  (2.27) 

where  (f)  is  the  particle-to-particle  angle  of  internal  friction,  and  a  is  the  cohesive 
force  (for  a  purely  frictional  contact  a  =  0).  In  the  terminology  of  crystallographic 
slip,  the  particle  contact  becomes  potentially  active  when  F  =  0.  Equation  (2.27) 
represents  the  yield  condition  corresponding  a  straight  line  separating  the  elastic 
and  plastic  responses.  The  elastic  region  is  given  by  the  set 

E={fT£^^  \F:=\fT\-a  +  fNt^ri(j><0}  .  (2.28) 


Following  plasticity  theory,  the  flow  rule  defines  the  direction  of  the  plastic  slip 
vector  by  assuming  a  plastic  potential  function  of  the  form 


G  =\fT\- a  +  /at  tan  V', 


(2.29) 


where  the  friction  angle  in  (2.28)  has  been  replaced  by  the  dilation  angle  Hence, 
the  flow  rule  G  in  which  =  0  gives 


jyp  = 


0 

1 


(2.30) 


Note  that  the  consistency  parameter  A  has  the  physical  meaning  of  being  the  particle 
contact  slip  itself. 


The  constitutive  equation  (2.26)  can  now  be  expanded  to  give 

In  terms  of  the  centroidal  particle  displacements  and  rotations,  the  indentation  rate 
8  for  small  8  can  be  found  from  (2.14)  as 


8  =  n^-{d^-  d^)  =  ■  {d^  -  d^) .  (2.32) 

Similarly,  the  slip  rate  7  for  small  8  is  found  from  (2.22)-(2.24)  to  be 

7  =  ,  (2.33) 

qAC  _-qA^  (2.34) 

^BC  _qB  ^  (2.35) 


in  which 


=  fea  •  X  — — —)  sec 9^  =  fes  •  x  — — —)  secO^  .  (2.36) 

V  +  tb  /  V  VA  +  rB  y 

The  hardening  law  represents  an  important  element  of  any  constitutive  relation. 
For  simplicity,  a  linear  hardening  law  has  been  assumed  of  the  form 

a  =  H'\jP\.  (2.37) 

Note  that  when  H'  =  0,  the  particle  contact  becomes  ideally  plastic.  Thus,  the 
consistency  condition  which  requires  that  the  yield  condition  be  satisfied  as  long  as 
the  contact  is  in  the  plastic  state  gives 

^=^'=('  +  f:) + 

The  loading/unloading  conditions  for  a  particle  contact  slip  system  can  now  be 
written  as  follows: 

i.  the  contact  unloads  or  is  inactive  (‘stick’  mode)  if 

F  <  0,  OT  F  =  0  and  A  <  0;  (2.39) 


ii.  the  contact  loads  or  is  active  (‘slip’  mode)  if 

jF  =  0  and  A  >  0. 


(2.40) 


The  constitutive  relations  (2.26)  can  be  summarized  by  relating  the  particle 
forces  F  to  the  contact  movement  V  as  follows: 


where 


and 


F  =  H{6)fC^PV , 


VP  =  v-{i-k:), 


1C  = 


H' 


(2.41) 


(2.42) 


(2.43) 


and  1  is  the  second-rank  identity  tensor  and  are  the  elastic-plastic  moduli.  The 
constitutive  relation  (2.41)  for  a  general  particle  contact  has  two  branches  which 
corresponds  to  continued  loading  (‘slip’)  and  to  unloading  (‘stick’)  (see  Figure  2). 
Clearly,  if  the  slip  system  is  inactive  or  unloads  (‘stick’),  the  response  is  purely 
elastic  and  (2.42)  reduces  to  =  fC^. 


2.4.  Contact  force  integration 

Let  the  superimposed  V  denote  the  time  rate  of  change  co-rotational  with  the 
contact  normal.  Consider  a  rate  constitutive  equation  of  the  form 

(2.44) 

where  R  represents  a  transformation  operator  which  rotates  the  local  contact 
element  axes  to  the  global  (Cartesian)  coordinate  system.  The  first  term  on  the 
right-hand  side  of  (2.44)  is  given  by  (2.41)  and  represents  the  material  response  due 
to  deformation  of  the  contact  element.  The  second  term  accounts  for  rotational 
effects.  For  two-dimensional  problems  on  the  plane  xi,X2,  the  rotation  and 
rotational  rate  take  the  form 

cos  0  —  sin  0 

sin  0  cos  0 


R  = 


(2.45) 


R=- 


sin  0  cos  0 
—  cos  0  sin  0 


0. 


(2.46) 


where  0  is  the  orientation  of  the  contact  normal  with  respect  to  the  positive  Cartesian 
a;i-axis.  Since  N  is  fixed  and  0^  represents  the  rotation  of  n  relative  to  N,  0  =  0^, 
see  (2.36). 

The  form  of  (2.44)  facilitates  an  exact  analytical  integral  of  forces  which  leads 
to  an  incremental  constitutive  equation  of  the  form 


A/ 


=  / 

Jin 


in+1  V 

^  dt  =  Rn+l^n+l 


Rn^n  —  fn-^\  fn  i 


(2.47) 


where  fn+i  —  Rn+l^ fn  —  Rn^ n.1  s-ud  ^ —  {f N 1  Jt}  •  The  foices  fn  and 
will  be  used  in  Chapter  4  to  construct  internal  force  vectors  for  each  contact 
element. 
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Expressions  for  the  contact  forces  ^n+l  will  be  obtained  by  considering 
contacting  particles  A  and  B.  Assume  that  at  time  the  positions  of  the 
particle  centroids,  x'^  and  ,  are  well  defined.  Then  integrating  (2.41)  yields 
the  component  /jy: 

(/7v)n+l  =  -He{Sn+l)kNSn+l  ;  ^n+1  =  +  r"®  -  ||in+l||  ,  (2.48) 

where 

j„+i  =  +  df+i  -  ) ,  (2.49) 

and  He  is  a  ramp  function  of  the  form 

(  1,  if  ^n+i  >  £; 

He{^n+\)  —  ■!  ^n+l/S)  if  0  <  <  £;  (2.50) 

I  0,  if  (5„+i  <  0. 

in  which  e  is  a  ‘sufficiently  small’  number.  This  ramp  function  approximates  the 
Heaviside  function  such  that  as  e  ^  0,  He{S)  — >  H{8).  It  makes  the  integrated 
constitutive  equation  continuous  in  the  sense  of  Lipschitz  [33]  and  regularizes  an 
otherwise  non-regular  contact  problem.  The  ramp  function  renders  the  discrete 
contact  forces  amenable  to  linearization. 

The  tangential  force  /y  reflects  the  previous  slip  mode  between  two  contacting 
particles  and  contains  the  material  memory.  Thus,  assuming  that  at  time  the 
magnitude,  and  the  sense  of  the  force  {fT)n  are  known,  the  slip  increment  can  then 
be  simply  evaluated  as 

A7  =  A  ,  (2.51) 

where 

A0^^  =  Ae^  -{6^+1- e^),  (2.52) 

Ad^^  =  A0^  (2.53) 

and 

AO^  =  sin“^(e3  •  x  =  sin"^(e3  •  nf  x  n^+i) .  (2.54) 

Now,  the  contact  modes  (2.39)  and  (2.40)  can  be  defined  for  a  finite  incremental 

motion  assuming  the  particles  remain  in  contact  during  the  time  step  as  follows; 
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i.  the  contact  is  in  ‘stick’  mode  if 


I  (Mn  +  I  <  Q;„  +  tan  (j)  , 

and 

/t  =  (/r)n  +  Hekr^'j , 

A7P  =  0 , 

otherwise; 

ii.  the  contact  is  in  ‘slip’  mode: 

fx  =  sign(/r  )an+i  -  (/jv)n+l  fan  ^ , 

where 

/t  “  (/r)n  +  Hskx^'y  , 

On+l  =  +  ^'|A7^|  , 

lA  PI  l/rl -Q^n  +  /;ytan<^ 

'  '  H,kT  +  H' 


(2.55) 


(2.56) 

(2.57) 


(2.58) 

(2.59) 

(2.60) 

(2.61) 


and  represents  the  trial  tangential  force,  On+i  denotes  the  updated  cohesive 
force,  and  A7P  is  the  plastic  slip. 

The  steps  necessary  to  calculate  /jv  and  /j  are  summarized  in  Box  1.  Note 
that  because  the  load  increment  is  finite  and  the  particle  assembly  may  be  locally 
discontinuous,  it  may  happen  that  during  the  iteration  a  new  contact  is  formed 
between  two  previously  disconnected  particles.  Unless  a  very  small  load  increment 
is  employed,  it  may  never  be  known  when  the  initial  contact  is  actually  formed 
between  these  two  particles.  Box  1  thus  presents  two  options  to  bracket  the  expected 
response.  For  lOPT  =  1,  initial  contact  is  assumed  to  take  place  at  time  fn+i,  in 
which  case,  {fT)n+l  becomes  identically  equal  to  zero  (see  Step  3);  for  lOPT  =  0, 
initial  contact  is  assumed  to  take  place  instantaneously  at  time  i„. 
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Box  1.  Contact  force  calculations. 


Step  1.  Compute  ^n+l,  {fN)n+l  =  -He{Sn+l)kNSn+l- 
Step  2.  If  (/jv)»+i  >  0,  (/iv)n+i  =  (/T)n+1  =  0  and  return. 

Step  3.  If  =  0  and  lOPT  =  1,  (/T)n+i  =  0  and  return. 

Step  4.  Compute  A7. 

Step  5.  Compute  (/T)n+i  =  (/r)n  +  HekxA'y. 

Step  6.  Compute  /max  =  On  -  (/jv)n+i  tan  /. 

Step  7.  If  |(/r)n+i|  >  /max,  (/r)n+i  =  “  (/iv)n+l  tan /)  sign(/r)„+i . 

Step  8.  Return.  □ 
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Chapter  3 

Micromechanical-Macromechanical  Connections 


3.1.  Introduction 

The  particulate  nature  of  granular  materials  causes  the  macroscopic  applied 
loads  to  be  carried  at  contact  points  between  grains.  Thus,  employing  continuum 
concepts  to  describe  the  overall  response  of  a  material  with  microstructure  and 
particulate  mechanics  to  describe  of  constituent  particle  motions  presents  a  modeling 
challenge.  For  example,  since  stress  represents  a  macroscopic  continuum  concept, 
its  use  for  a  particulate  medium  requires  careful  micromechanical  considerations  of 
particle  contact  force  transmission.  Likewise,  the  notion  of  macroscopic  kinematic 
quantities  demands  special  interpretation  of  the  role  micromechanics  plays  in 
describing  the  motion  of  individual  particles.  Hence,  the  overall  stress  and  kinematic 
relationships  establish  crucial  connections  between  the  macromechanical  description 
and  the  fundamental  underlying  particulate  mechanics. 

This  chapter  describes  the  connection  between  the  particle  contact  forces  and 
the  overall  stress.  The  kinematical  description  of  the  overall  velocity  gradient  will 
also  be  considered  in  terms  of  an  average  of  appropriate  microscopic  quantities. 
These  two  essential  connections  then  form  the  basis  for  the  macroscopic  constitutive 
relation.  The  development  of  the  criterion  for  predicting  the  onset  of  localized 
deformation  follows  naturally  from  this  relation. 


3.2.  Macroscopic  stress  fieLd 

The  overall  stress  acting  on  a  particle  assembly  produces  contact  forces  at 
contacting  points  of  the  constituent  particles.  Independent  of  the  nature  of  these 
contacts,  the  basic  objective  in  this  section  is  to  describe  the  overall  stress  in  terms 
of  the  contact  forces  and  some  geometric  characteristics  of  the  granules. 

Consider  a  representative  sample  of  a  granular  mass  of  volume  V  and  surface 
area  S.  Let  two  control  volume  particles  A  and  B,  with  centroids  at  and 
respectively,  have  a  contact  point  at  x^^  as  shown  in  Figure  3.  The  branch  vector 
connects  the  centroid  of  particle  A  to  the  centroid  of  particle  B  and  is  given  by 

l^B  =  x^-xA.  (3.1) 

Let  fA^  and  f^A  denote,  respectively,  the  contact  forces  exerted  on  grain  A  by 
grain  B  and  vice  versa,  and  therefore  fA^  =  —f^A 

Neglecting  inertia  and  gravity  terms  the  balance  of  linear  momentum  for 
particle  A  requires 

=  (3.2) 

/)=] 

where  fA^  represents  the  contact  forces  exerted  on  particle  A  by  particle  and 
K  denotes  the  so-called  coordination  number  for  particle  A  (i.e.,  the  number  of 
particles  contacting  A). 

The  balance  of  angular  momentum  for  particle  A  yields 

^  fAP  X  {xAP  -xA)  =  0,  (3.3) 

where  is  the  position  vector  of  the  contact  point  between  particles  A  and 
Summing  the  contributions  from  all  particles  contained  in  the  control  volume  gives 

N 

^/«xr=0,  (3.4) 

^=1 
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Figure  3.  Two  control  volume  particles  with  centroids  at 
and  contact  point  at  x^^,  and  branch  vector 

where  a  represents  a  particle  contact  such  as  AB^  and  N  denotes  the  total  number 
of  contacts  in  the  control  volume.  An  important  implication  of  equation  (3.4)  is 
that  it  renders  the  following  tensor  symmetric 

N  N 

=  (3.5) 

0  =  1  0  =  1 

To  relate  the  contact  forces  to  the  macroscopic  stress,  the  principle  of  virtual 
work  is  employed  in  the  manner  of  Christoffersen  et  al.  [29].  First,  the  overall  stress, 
a  is  given  by  its  natural  definition  [24,  25] 

a-  =  ^  a{x)  dV ,  (3.6) 

where  (^{x)  represents  the  variable  stress  field  in  equilibrium  with  the  overall  applied 
boundary  traction  T,  i.e. 

a-  •n  =  T  on  5 ,  (3.7) 

where  h  denotes  the  exterior  unit  normal  to  S. 

Now,  consider  a  sufficiently  smooth  overall  virtual  displacement  u  that  results 
in  a  virtual  separation  of  the  ath  contact  force.  Assuming  that  the  traction  T 
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is  in  equilibrium  with  the  contact  forces  the  virtual  work  principle  requires 


]_ 

V 


r  ^ 

/  T  •  w  d5  =  ^ 

a=l 


(3.8) 


Since  the  virtual  displacement  field  can  be  chosen  arbitrarily,  consider  the  linear 
field 


u  =  <f>  ■  X  +  c, 


(3.9) 


where  <j}  and  c  denote,  respectively,  an  arbitrary  second  order  tensor  and  an  arbitrary 
vector.  The  virtual  contact  separation  compatible  with  this  field  is  given  by 


.  r . 


(3.10) 


Appendix  A  shows  that  this  representation  of  (f)  ■  l°‘  is  reasonable  to  a  first  order 
approximation  to 


Substituting  (3.9)  and  (3.10)  into  (3.8)  and  employing  the  divergence  theorem 


gives 


cf>: 


N 

dF  -  ^  (8)  r 

a=l 


=  0. 


(3.11) 


With  the  symmetry  property  (3.5)  and  the  definition  of  the  overall  stress,  equation 
(3.11)  reduces  to 

1  ^ 

&  =  5l](r  ®i“+i“®/").  (3.12) 

a=l 

Note  that  this  expression  is  dimensionally  consistent  if  N  is  interpreted  as  the 
number  of  contacts  per  unit  volume.  The  derivation  of  the  macroscopic  stress  & 
requires  the  assumption  that  the  virtual  displacements  are  small.  Thus,  the  overall 
stress  inherits  a  spatial  definition,  and  &  may  be  termed  the  overall  Cauchy  stress 
arising  from  the  contact  forces. 


It  is  important  to  recognize  that  from  each  pair  of  contact  forces  at  a  given 
contact  (i.e.,  or  for  contact  AB)  only  one  enters  the  summation  in  (3.12). 
The  choice  of  the  vector  determines  the  choice  of  (i.e.,  when  I®  =  , 

jnanner,  at  each  contact,  a  pair  of  vectors  Z®  and  /®  has  a 
unique  tensor  product  and  contribution  to  the  overall  stress. 


3.3.  Macroscopic  deformation  field 


In  this  section,  the  relationships  between  the  macromechanical  and 
micromechanical  kinematic  quantities  will  be  described.  Consider  the  same 
representative  control  volume  composed  of  discrete  particles,  and  again  let  V  denote 
the  volume  and  S  the  surface  area.  Recall  that  this  volume  represents  a  point  in 
the  macroscopic  sense  and  that  the  overall  quantities  (stresses,  deformations,  etc.) 
are  the  macroscopic  responses  themselves.  Now,  assume  a  velocity  field  v  E 
that  is  sufficiently  smooth  over  the  problem  domain,  i.e.. 


dv 

L  =  -7—  =  constant  in  V , 
ox 


(3.13) 


where  V  =  V  U  S'  is  the  closure  of  V,  and  L  is  the  uniform  overall  velocity  gradient. 
Note  from  the  overall  velocity  gradient  Z/,  the  overall  rate  of  deformation  D  and 
spin  W  can  be  evaluated  via 

£)  =  symm(L);  W  =  skew(Z) ;  L  =  D  +  W.  (3-14) 

A  micromechanical  connection  with  (3.13)  and  (3.14)  can  be  described  in  the  manner 
of  Nemat-Nasser  and  Mehrabadi  [23]. 

Consider  a  unit  cell  centered  at  a  point  with  position  vector  x  experiencing  a 
uniform  velocity  gradient  L.  As  this  unit  cell  shrinks  to  a  point,  one  recovers  the 
macroscopic  definition  (3.13).  However,  the  cell  must  be  large  enough  to  contain  a 
sufficient  number  of  particles  to  capture  the  overall  behavior  of  the  material.  Within 
the  control  volume  the  local  velocity  field  v  —  v(x)  will  not  generally  be  a  linear 
function  of  x  but  will  vary  according  to  the  typically  irregular  particle  motions.  As 
a  result,  the  associated  velocity  gradient  L  will  also  be  irregular  and  is  given  by  the 
expression 


L  = 


dv 


(3.15) 


where  x  represents  the  position  vector  of  a  point  in  the  neighborhood  of  x  within 
the  unit  cell. 


Applying  the  divergence  theorem  to  (3.13)  and  (3.15)  yields 


(3.16) 
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and 


where  the  symbol  (•)  denotes  an  average  over  the  control  volume,  and  n  again 
represents  the  exterior  unit  normal  vector  to  S. 

If  the  boundary  particles  move  uniformly  such  that  v  =  u  on  S',  then  the 
left-hand  sides  of  (3.16)  and  (3.17)  are  identical,  which  implies  that 

I=(i).  (3.18) 

The  physical  meaning  of  this  equation  is  that  the  macroscopic  velocity  gradient 
equals  the  volume  average  of  the  local  particle  velocity  gradients. 

Now,  let  Ait  be  a  sufficiently  smooth  macroscopic  displacement  field  associated 
with  the  overall  velocity  field  i),  measured  with  respect  to  the  configuration  at  time 
station  Then,  the  updated  position  vector  at  any  time  t  is 

x  =  Xn  +  Au,  (3.19) 


where  is  the  configuration  at  the  beginning  of  the  time  step. 

The  associated  overall  deformation  gradient  for  this  motion  is 

-  dx  dAu 

F  =  —  =  Fn  + 


dX 


dX  ’ 


(3.20) 


where  Fn  =  dXnjdX.  Now,  if  is  taken  as  the  reference  configuration,  then 
(3.20)  degenerates  to 


dx  ^  dAu  -  -1 

^  dXn  dXn~ 

Since  the  overall  velocity  gradient  L  is  a  spatial  tensor,  then 

L  =  F 

^hFn-K'  r' 

=  hf~\ 
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(3.21) 


(3.22) 

(3.23) 

(3.24) 


i.e.,  L  does  not  depend  on  the  choice  of  the  reference  configuration. 

Equation  (3.19)  can  be  used  to  obtain  the  updated  configuration  at  time  tn+i 
as 

-X^n+l  =  d"  •  (3.25) 

Now,  if  Auji  is  assumed  to  be  a  homogeneous  function  of  degree  one,  then  Euler’s 
theorem  for  homogeneous  functions  gives 

=  A  •  ,  (3.26) 

where  A  =  dAunIdXn-  The  macroscopic  velocity  field  associated  with  (3.26)  is 

v  =  A- Xn.  (3.27) 

Substituting  (3.26)  and  (3.21)  into  (3.24)  gives  the  overall  velocity  gradient  at  time 
as 

jb  =  A  •  (1  +  A)-^ .  (3.28) 

Hence,  for  L  to  have  a  uniform  distribution  over  the  control  volume  V,  the  tensor  A 
should  be  spatially  constant  over  V.  Physically,  the  macroscopic  tensor  A  provides  a 
measure  of  the  overall  deformation  experienced  by  the  macroscopic  point  of  interest. 

Conceptually  equation  (3.26)  indicates  that  if  adjacent  control  volumes  lie 
sufficiently  close  to  each  other,  they  may  overlap.  This  observation  is  consistent 
in  order  for  the  macroscopic  tensor  A  (and  thus,  the  overall  velocity  gradient  L)  to 
have  a  sufficiently  smooth  distribution  over  the  problem  domain. 


3.4.  Macroscopic  constitutive  relation 

A  macroscopic  constitutive  relation  will  be  useful  for  characterizing  the  overall 
behavior  of  the  granular  material  in  terms  of  its  particulate  nature.  As  will  be 
seen  in  the  next  section,  the  macroscopic  constitutive  relation  plays  a  central  role 
in  the  development  of  the  localization  criterion  for  the  material.  Localization  is 
interpreted  herein  as  an  instability  in  the  macroscopic  constitutive  description  of 


the  inelastic  deformation  of  a  material.  The  constitutive  equations  will  also  be  of 
particular  importance  in  the  development  of  the  response  of  a  macroscopic  point  to 
an  applied  overall  stress  history  presented  in  Chapter  5. 

The  macroscopic  constitutive  relation  depends  heavily  upon  the  micromechan- 
icahmacromechanical  connections  that  bridge  the  continuum  and  particulate 
descriptions  of  the  granular  material.  First,  the  dependency  of  the  overall 
stress  on  these  connections  will  be  determined.  The  particle  contact  force 
integration  presented  in  Section  2.4  provides  the  contact  forces  for  evaluating 
the  macroscopic  stress  a  given  in  (3.12).  The  local  constitutive  relation  describes 
these  forces  in  terms  of  the  particle’s  centroidal  motion  which  in  turn  depends 
in  part  on  the  applied  overall  displacement  gradient  A.  In  terms  of  overall 
quantities,  an  implicit  macroscopic  constitutive  equation  can  be  written  to  capture 
this  relationship  as  follows: 

o-„+i  =  (3.29) 

i.e.,  the  overall  stress  is  a  function  of  the  components  of  the  displacement  gradient  A. 


In  principle,  (3.29)  can  always  be  written  in  rate  form, 
differentiating  (3.29)  with  respect  to  time  and  using  (3.28)  gives 

*-s- 

=  K  :L. 


For  example, 


(3.30) 

(3.31) 

(3.32) 


The  fourth  order  tensor  K  represents  an  overall  moduli  tensor  with  components 


Kijkl  =  ^  A  i^lm  +  ^Im)  i  (3.33) 

where  6im  denotes  the  Kronecker  delta.  When  the  components  of  the  displacement 
gradient  A  are  small,  then  equation  (3.33)  reduces  to 


dAki 


(3.34) 


Since  slip  governs  the  microstructural  rearrangement  mechanism,  a  vertex  will 
develop  on  the  macroscopic  yield  surface  [30,31].  In  other  words,  the  relation  (3.32) 


will  have  an  infinite  number  of  branches  corresponding  to  different  directions  of  L. 
Thus,  plastic  ‘normality’  in  conjugate  deformation  variables  is  lost  [30],  and  the 
moduli  tensor  K  is  restricted  to  only  minor  symmetry 

Kijki  =  Kjilk-  (3.35) 

The  minor  symmetry  of  K  with  respect  to  the  indices  ij  follows  directly  from  the 
symmetry  of  a.  The  minor  symmetry  with  respect  to  the  indices  kl  results  from  the 
fact  that  L  is  rotation-free,  i.e.,  from  a  macroscopic  standpoint  the  strains  are  small. 
Thus,  the  tensor  K  can  be  used  for  small-strain  stability  analysis  as  demonstrated 
in  the  next  section. 


3.5.  Localization  criterion 

Critical  conditions  are  now  sought  at  which  the  material’s  macroscopic 
constitutive  relation  may  allow  a  bifurcation  from  a  homogeneous  or  smoothly 
varying  deformation  into  a  highly  concentrated  non-uniform  deformation  in  a  planar 
band.  Outside  the  localized  band,  the  conditions  of  continuing  equilibrium  and 
continuing  homogeneous  deformation  must  be  met.  The  onset  of  localization  results 
from  the  loss  of  ellipticity  of  the  continuing  velocity  equilibrium  equation 

|i  =  0.  (3.36) 

The  problem  is  then  to  determine  conditions  at  which  (3.36)  is  satisfied  across  some 
planes  of  orientation  n. 

Consider  a  macros copically  homogeneous,  homogeneously  deformed  material 
subjected  to  quasi-static  increments  of  deformation.  For  the  velocity  field  to  remain 
continuous  at  bifurcation,  compatibility  requires  the  velocity  field  to  be  expressed 
as 

(3.37) 

where  superscripts  b  and  o  denote  band  and  outside  the  band  respectively,  and  g  is 
a  function  only  of  the  distance  across  the  band,  n  • »,  and  is  zero  outside  the  band. 
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The  condition  of  continuing  equilibrium  (3.36)  also  demands  that  the  traction  be 
continuous  across  the  discontinuity  planes 


n-  (T*  ■=  n-  <T  .  (3.38) 

With  (3.32),  the  constitutive  equation  inside  the  band  is 

a  =K^:i\  (3.39) 

and  substituting  in  (3.37)  and  (3.38)  yields 

(n-K^  ■nyg  =  n- (k^ -K^y.L\  (3.40) 

A  continuous  constitutive  response  is  assumed  prior  to  localization,  i.e.,  IT®  =  K^. 
Thus,  for  a  non-trivial  solution  for  g,  the  condition  for  localization  becomes 

det  [>4(n)]  =  det  [n  •  IT  •  n]  =  0 ,  (3.41) 


where  A  =  n  •  K  •  n  \s  the  so-called  ‘acoustic  tensor.’  Appendix  B  contains  the 
details  for  calculating  the  orientation  of  the  band  n  in  two  dimensions. 


Chapter  4 

Strain-driven  Problem 


4.1.  Introduction 

With  the  background  of  Chapter  2  and  Chapter  3,  a  solution  strategy 
may  be  constructed  to  derive  the  response  of  a  macroscopic  point  in  a  granular 
material  from  the  behavior  of  a  particle  assembly.  The  solution  model  adopts 
the  conventional  hypothesis  in  computational  plasticity  by  assuming  a  deformation 
driven  problem,  i.e.,  the  given  overall  displacement  gradient  determines  the  overall 
stress  response.  The  mathematical  formulation  prescribes  an  overall  uniform  motion 
to  an  assembly  of  particles  and  then  allows  the  particles  to  move  in  a  microscopic 
sense.  Throughout  the  particle  motions,  contacts  may  form  or  break  altering  the 
material’s  microstructure.  The  contact  forces  that  develop  between  particles  can 
then  be  used  to  evaluate  the  overall  stresses. 

This  chapter  considers  the  solution  of  the  deformation  driven  problem  in  detail. 
The  formulation  is  presented  in  the  context  of  an  analogy  with  the  finite  element 
method  where  the  control  volume  particles  represent  nodes  and  their  contacts  denote 
elements.  Special  attention  is  given  to  the  application  of  the  uniform  deformation 
to  the  assembly  of  particles.  The  notion  of  a  repeating,  or  periodic,  cell  will  be 
presented  which  allows  the  prediction  of  the  overall  response  of  the  control  volume 
without  the  imposition  of  boundary  particle  displacements. 


4.2.  Algorithm  for  problems  with  an  imposed  deformation  history 

Consider  a  two-dimensional  assembly  of  rigid  circular  particles.  Each  particle 
contains  two  translation  degrees  of  freedom  in  the  Xi,  X2  plane  and  one  rotational 
degree  of  freedom  about  the  xs  axis  as  shown  in  Figure  4.  The  constituent  particles 
form  a  moving  mesh  where  the  particles  have  been  replaced  by  nodes  and  the 
contacts  are  represented  by  the  stick/slip  elements  of  Chapter  2.  The  associated 
grid  continually  changes  as  particles  move  to  new  positions  breaking  contacts  and 
forming  new  ones,  i.e.,  severing  contact  elements  and  forming  new  ones. 

The  internal  force  vector  at  time  tn+i  for  a  typical  contact  element  number  ‘e’ 
connecting  particles  A  and  B  is  given  by 

/('exl)  =  ,  (4.1) 

where  the  balance  of  linear  momentum  (3.2)  requires  that  =  —fn+i-  Recall 

that  and  represent  that  internal  force  vectors  reckoned  with  respect  to 

the  Cartesian  axes  xi,  X2  and  the  rotation  matrix  R  relates  them  to  the  normal  and 
tangential  contact  forces  in  T  via 

/„+l  =  .?■»+. .  (4.2) 

where 

=  {/aTi/t}*  •  (4-3) 

The  corresponding  particle  degrees  of  freedom  are 

‘'(6K1)  =  ,  (4.4) 

where  represents  the  vector  of  unknown  nodal  displacements  for  particle  A 
and  denotes  the  unknown  rotation  of  particle  A  (see  Figure  4). 

Imposing  the  momentum  balance  equation  for  every  particle  in  the  control 
volume  and  neglecting  inertia  and  gravity  terms  results  in  an  equilibrium  equation 
of  the  form 

net 

FiNT{<L+l)=[Jr  =  0,  (4.5) 

e=l 
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Figure  4.  Particle  degrees  of  freedom  on  the  Xi ,  X2  plane. 


where  the  symbol  U  denotes  the  assembly  operator.  To  determine  the  vector  of 
nodal  displacements  that  satisfies  (4.5)  at  each  time  step,  Newton’s  method 

with  line  search  can  be  employed: 

f’W(<+l)A<i=F„r(<+i);  <++■,- <+i-«Ad;  <?,=<<.,  (4.6) 

where  0  <  a  <  1  is  a  line  search  parameter  introduced  to  insure  that  the  iteration 
(4.6)  is  norm-reducing  [34,35].  In  (4.6),  the  algorithmic  tangent 

operator  derived  from  the  element  contribution 

net 

=  [J  ■  (4-7) 

e=l 

By  definition,  the  solution  has  converged  when 

ll-^/Vrll/ll-^/Vrll  <  (4-8) 


where  rtol  represents  a  prescribed  error  tolerance. 


Taking  the  derivatives  of  (4.1)  results  in  the  following  element  contributions 
to  the  tangent  operator  (omitting  subscripts  ‘n  +  1’  for  time  step  and 
superscripts  ‘i’  for  iteration  counter,  for  simplicity): 


no'') 

f(0B)  - 

r^Srid^) 

-f(d^) 

-f'(d^) 

-/(««) 

r^Srid'^) 

r»/i,(e4) 

r‘>ST(d‘‘) 

(4.9) 
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Box  2.  Strain- driven  algorithm. 


i.  Given  An,  and  i  =  0. 

ii.  Compute 

iii-  ~  ^n+l  “  [-^INTi^h+l)]  ’’*• 

iv.  r»+i  =  U:i'ir(<5-\)- 

V.  IF  ||r*+i||/||r0l|  <rtol,  GO  TO  5. 
vi.  i  ?  -f  1  and  GO  TO  ii.  □ 


where 


f(d^)  =  -f'{d^) ;  r-8/(«-<)  =  ;  (4.10) 

flpid'')  =  -f!r{d‘>) ;  r‘>flr{0^)  =  .  (4.11) 

Appendix  C  contains  a  complete  discussion  of  the  derivatives  contained  in  (4.9). 

In  general,  the  matrix  Fif^Ti^n+i)  neither  symmetric  nor  positive 

definite.  Furthermore,  it  may  be  singular  with  isolated  particles  or  clusters  of 
particles.  The  matrix  na^y  have  a  changing  bandwidth  from  iteration 

to  iteration  as  a  consequence  of  new  contacts  being  formed  and  old  contacts  being 
broken.  Section  4.4  discusses  these  numerical  issues  in  detail. 

Box  2  contains  the  steps  necessary  for  solving  a  deformation-driven  problem. 
The  algorithm  depends  heavily  upon  the  calculation  of  the  internal  force  vector 
f^{dn-\-i)  which  represents  the  particle  contact  force  contribution  to  the  internal 
force.  Recall  that  Chapter  2  described  the  role  that  the  spring  stiffnesses  kjif  and 
kx  play  in  the  contact  constitutive  equations  used  for  calculating  the  contact  forces. 
In  the  context  of  the  deformation-driven  problem  of  Box  2,  these  spring  stiffnesses 
may  be  considered  as  penalty  parameters  rather  than  as  physical  measures  of  the 
true  particle  rigidities. 


4.3.  The  notion  of  a  repeating  cell 


The  dependence  of  the  vector  on  the  imposed  displacement  gradients  A 
follows  directly  from  the  initial  imposition  of  a  uniform  motion  of  the  particle 
centroids  according  to  (3.25)  and  (3.26).  Assume  that  the  particle  centroids  are 
given  by  =  1,2,  N,  then  the  new  configuration  of  a  typical  particle  A 

can  be  written  as 

®^  =  (l  +  A)•X^  (4.12) 

This  motion  will  generally  perturb  the  momentum  balance  equations  as  contacting 
particles  either  separate,  overlap,  or  slip  during  the  uniform  motion.  The  unbalanced 
forces  create  residuals  in  the  force  vector, 

r  =  FiNT{A-X^)^0,  (4.13) 

where  (4.12)  has  been  used  as  an  initial  estimate  for  Newton’s  method  (4.6).  The 
continued  dependence  of  on  A  follows  from  the  dissipation  of  the  residual 
vector  r.  Through  the  Newton  iteration,  the  particles  seek  their  equilibrium 
configuration  while  the  sides  of  the  control  volume  are  held  fixed  according  to  the 
displacements  produced  by  A.  Thus,  the  sides  of  the  unit  cell  constrain  the  particles 
to  move  according  to  the  deformed  configuration  imposed  by  A. 

As  the  iterations  progress,  the  particles  will  translate  and  rotate  to  find  their 
equilibrium  positions,  which  could  be  drastically  different  from  the  imposed  uniform 
motion.  Hence,  it  may  happen  that  v  ^  v  even  on  S,  and  so  one  cannot  simply 
prescribe  the  motion  of  the  boundary  particles  in  order  to  constrain  the  problem.  To 
determine  the  particle  motion  that  is  independent  of  the  boundary  displacements, 
the  notion  of  a  repeating,  or  periodic,  cell  has  been  employed. 

Periodicity  of  the  unit  cell  requires  that  the  control  volume  of  interest  be 
surrounded  in  all  directions  by  identical  parallelepipeds  (or  rhombuses,  in  two 
dimensions).  Physically,  the  notion  of  a  periodic  cell  means  that  if  a  particle  leaves 
the  control  volume  an  identical  mirror  image  particle  will  enter  the  unit  cell  from 
a  contiguous  unit  cell.  Periodicity  can  be  enforced  for  each  pair  of  potentially 
contacting  particles,  A  and  B,  by  assuming  that  particle  A  could  be  in  contact  with 


Figure  5.  Periodic  control  volume  subjected  to  a  uniform 
deformation.  As  a  particle  A  leaves  the  reference  cell  an 
identical  mirror  image  particle  A'  enters  the  ceU. 


either  B  or  one  of  the  images  of  jB  in  a  neighboring  parallelepiped  (or  rhombus). 
Figure  5  shows  a  graphical  representation  of  how  periodicity  may  be  imposed  for 
implicit  two-dimensional  plane-strain  calculations. 

In  order  to  describe  the  basic  repeating  unit  cell,  let  the  basis  vectors 
represent  the  sides  of  a  parallelepiped  (or  rhombus)  that  defines 
the  control  volume  V  in  the  reference  configuration.  Assume  that  one  corner  of 
the  parallelepiped  (or  rhombus)  coincides  with  the  origin  of  the  reference  frame  so 
that  the  ^^’s  become  the  position  vectors  of  the  corners  of  V .  In  addition,  let  the 
vectors  (yjj,  ^^25  •  •  •  5  V’njd)  represent  the  sides  of  the  same  volume  in  the  current 
configuration  such  that 

(Pi  =  {1  + A) i  =  l,2,...,ngd-  (4.14) 

Now,  let  A  and  B  represent  two  potentially  contacting  particles  with  respective 
centroids  at  and  measured  with  respect  to  the  current  configuration.  With 
respect  to  the  baisis  vectors  tpi,  the  branch  vector  connecting  the  particle  centroids 
x'^  and  x^  can  be  written  as 

x^  —x^  =  ai(f>i ,  (4.15) 
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Repeating 
/  Cell 


Figure  6.  Enforcement  of  the  periodicity  condition:  contacts 
are  checked  between  A  and  B,  and  between  particle  A  and 
the  three  image  particles  B'. 

where  a  sum  is  implied  on  i  =  1,2, ...  ,nsd,  and  the  q;,’s  are  some  determinate 
multipliers.  Since  the  yj,’s  form  a  basis,  then  the  Oi’s  are  all  distinct. 

Now,  periodicity  of  the  cell  can  be  imposed  by  assuming  that  particle  B  and  its 
images  are  located  at  the  corners  of  a  parallelepiped  (or  rhombus)  that  is  identical  to 
the  unit  cell  of  interest  and  inscribes  the  centroid  of  particle  A  as  shown  in  Figure  6. 
In  terms  of  the  current  configuration  vectors  <pi,  the  centroidal  coordinates  of  particle 
B  and  its  images  are  given  by 

x^  =  x^-^i(pi,  (4.16) 

where  a  sum  is  again  implied  on  i,  and  the  ^i’s  are  permutations  of  (0,  sign(Q!,)).  For 
example,  when  l3i  =  02  =  •••  =  0n,i  =  0,  then  x  =  x^ ,  which  is  the  centroidal 
position  vector  for  particle  B',  when  0i  =  sign(ai)  and  02  =  ■■■  =  0n,ci  —  0,  then 
X  —  —  sign(ai)(^i,  which  is  the  centroidal  position  vector  for  an  image  of  B, 

and  so  on.  Figure  6  shows  a  graphical  representation  of  the  enforcement  of  the 
periodicity  condition  for  implicit  two-dimensional  plane-strain  calculations. 

Note  that  if  particle  A  lies  in  contact  with  B,  it  can  never  be  in  contact  with  any 
of  the  images  of  B.  The  element  contributions  to  the  residual  force  vector  r  represent 


particle  contact  element  contributions  and  not  particle  contributions.  Thus,  one 
and  only  one  contact  force  vector  can  be  added  to  the  residual  regardless  of  the 
multiplicity  of  the  particle  images. 

An  important  feature  of  the  periodic  cell  lies  in  its  ability  to  allow  the  prediction 
of  the  macromechanical  response  without  the  imposition  of  boundary  particle 
displacements.  However,  since  periodic  cells  inherently  contain  zero-energy  modes, 
the  periodicity  condition  alone  does  not  guarantee  a  fully  constrained  boundary- 
value  problem.  Rigid-body  modes  may  be  eliminated  by  fixing  the  motion  of  any  one 
particle  in  V,  provided  that  this  particle  belongs  to  the  ‘principal  force  chain.’  Rigid- 
body  modes  cannot  be  removed  by  fixing  an  isolated  particle,  nor  any  particle  in  an 
isolated  cluster.  The  next  section  discusses  numerical  examples  and  implications  of 
the  periodic  cell  and  zero-energy  modes. 


4.4.  Numerical  simulations 

This  section  presents  the  results  of  two-dimensional  plane-strain  simulations  on 
granular  assemblies  composed  of  either  regular  or  random  initial  packings  of  circular 
disks.  Some  fundamental  properties  of  granular  materials  are  demonstrated  such  as 
softening  and  anisotropy.  Numerical  difficulties  encountered  during  the  simulation 
process  are  reported  as  well.  Unless  otherwise  stated.  Table  1  defines  the  model 
parameters  used  in  all  analyses.  The  overall  stresses  have  been  calculated  using  a 
volume  equal  to  the  initial  cross  sectional  area  Aq  times  a  unit  thickness.  The  error 
tolerance  rtol  used  in  (4.8)  ensures  that  the  iteration  has  converged  sufficiently,  thus 
minimizing  the  propagation  of  numerical  errors. 

Sixty  four-particle  regular  assembly 

The  initial  granular  assembly  configuration  is  shown  in  Figure  7.  fn  this 
example,  the  assembly  is  assumed  to  be  composed  of  64  uniform  circular  particles  of 
radius  r  =  1.0  units  and  arranged  initially  in  simple  cubic  packing  except  that  the 
adjacent  rows  of  particles  have  been  shifted  by  5  degrees.  The  particle  on  the  lower 


Table  1.  Model  parameters. 


Error  tolerance:  rtol  =  1.0  x  10“®. 

Normal  spring  stiffness:  kj\f  =  1.0  x  10^. 
Tangential  spring  stiffness:  kx  =  1.0  x  10^. 
Particle  friction  angle:  (f)  —  30°. 

Particle  contact  cohesion:  ao  =  0. 

Contact  hardening  parameter:  H'  =  100. 
Ramp  function  parameter:  e  =  0.10. 


left-hand  corner  was  fixed  against  translation  and  rotation  to  arrest  the  zero-energy 
modes  present  in  the  assembly. 

The  assembly  was  then  subjected  to  combined  isotropic  compression  and  shear 
by  prescribing  ten  increments  of  motion,  each  increment  defined  by  the  following 
elements  of  the  tensor  operator  A  in  (3.26):  An  =  A22  =  —1.0  x  10“^, 

Ai2  =  A21  =  —1.0  X  10“^.  If  second-order  deformations  are  ignored,  then 

the  total  motion  corresponds  to  normal  compressions  of  en  =  €22  =  —1%  and 
tensorial  shear  strains  of  ei2  =  £21  =  —10%.  The  final  deformed  configuration 
after  running  the  analysis  over  ten  time  steps  appears  in  Figure  8. 

Figure  8  shows  that  during  the  course  of  deformation,  rows  of  particles  have 
separated  into  four  isolated  clusters.  This  “strain-softening”  effect  is  a  typical 
result  when  a  granular  medium  with  an  initially  collapsible  structure  is  subjected 
to  shearing  deformation  which  is  far  beyond  the  material’s  ability  to  compact. 
Numerically,  this  phenomenon  will  manifest  in  the  form  of  zeros  appearing  on  the 
diagonal  of  the  factorized  tangent  operator,  which  are  not  operated  upon  during  the 
back  substitution  process. 

A  different  result  can  be  obtained  with  a  slightly  altered  strain  history.  Using 
the  same  initial  configuration  of  Figure  7,  the  following  overall  uniform  motions  are 


now  prescribed;  one  time  step  of  An  =  A22  =  — 1-0  x  10“^,  An  =  A21  =  0; 
followed  by  one-hundred  time  steps  of  A21  =  — 1.0  x  10~®,  An  =  A22  =  A12  =  0. 
If  second-order  effects  are  ignored,  this  motion  corresponds  to  a  total  volumetric 
strain  of  cn  -f  €22  =  —2%  followed  by  a  sidesway  of  10%  (which  produces  simple 
shearing  and  rigid  body  rotation).  Results  of  the  simulations  are  shown  in  Figure  9 
and  Figure  10. 

Figure  9  shows  the  final  deformed  configuration  characterized  by  a  stable 
microstructure.  The  isotropic  compression  produces  prestressing  effects  on  the 
elastic  springs  and  prevents  the  particle  contacts  from  breaking  during  the  shearing 
process.  Since  shearing  involved  essentially  horizontal  particle  translation  and 
scraping,  the  overall  response  of  the  assembly  is  a  direct  function  of  the  tangential 
spring  constant  kx-  Thus,  for  a  constant  H'lkx,  the  overall  stresses  can  be 
normalized  with  respect  to  the  spring  constant  kx- 

In  Figure  10,  the  normalized  overall  shear  stress,  012! kxi  is  plotted  versus  the 
overall  tensorial  shear  strain,  eu-  Specifically,  the  overall  shear  stress  is  evaluated 
from  (3.12)  as 

^12 = +  'r/f>  (^-17) 

The  initial  straight-line  portion  of  the  overall  stress-strain  plot  represents  the  elastic 
stretching  of  the  tangential  springs,  and  is  thus  also  normalizable  with  respect  to  the 
spring  constant  kx-  On  the  other  hand,  decreasing  the  moduli  ratio  H' /  kx  decreases 
the  cri2//i:T'response  at  post-yield.  The  5-degree  initial  offset  between  adjacent  rows 
of  particles  causes  sequential  yielding  at  contact  points,  thereby  resulting  in  non- 
uniform  tangential  shear  moduli  at  post-yield. 

Sixty-particle  periodic  assembly 

The  initial  configuration  for  this  example  is  the  same  as  in  Figure  7  but  with 
four  interior  particles  removed  so  that  the  resulting  assembly  is  represented  by  four 
repeating  cells.  The  entire  assembly  is  subjected  to  50  increments  of  motion,  each 
increment  defined  by  the  following  elements  of  A:  An  =  A22  =  —4.0  x  10“^, 
A21  =  Ai2  =  1-0  X  10~^.  If  second-order  effects  are  ignored,  then  the  total  motion 
is  equivalent  to  a  total  volumetric  strain  of  cn  -b  €22  =  —2%  and  a  total  tensorial 
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shear  strain  of  ei2  =  e2i  =  5%.  The  resulting  deformed  configuration,  shown  in 
Figure  11,  shows  the  expected  periodicity  exhibited  by  the  deformed  assembly.  If 
this  problem  had  been  analyzed  using  only  one  of  the  four  repeating  cells,  then  the 
macroscopic  stresses  obtained  from  the  sum  of  the  contact  forces  for  each  reduced 
cell  would  have  been  one-fourth  of  those  obtained  from  the  full  unit  cell. 

Fifty  eight-particle  irregular  assembly 

The  initial  configuration  of  the  control  volume  for  this  example  is  shown  in 
Figure  12.  Here,  the  repeating  cell  resembles  the  one  used  in  the  sixty  four-particle 
simulation  except  that  six  of  the  original  particles  in  the  assembly  were  removed. 
The  control  volume  was  subjected  to  38  incremental  motions,  each  increment 
described  by  the  following  elements  of  A\  An  =  An  =  — 7  x  10~^  and 

An  —  An  =  1  x  10“^.  After  the  final  load  step,  the  deformed  configuration  is 
shown  in  Figure  13. 

Figure  13  shows  the  deformed  configuration  for  the  case  where  no  particle  in 
the  assembly  was  fixed  against  rigid-body  motion.  When  this  problem  was  run  with 
the  particle  on  the  lower  left-hand  corner  fixed  against  translation  and  rotation,  then 
the  results  were  the  same  except  for  the  rigid-body  modes.  Thus,  provided  that  the 
assembly  remains  stable,  rigid-body  modes  may  be  allowed  without  afflicting  the 
computed  overall  stress-strain  response.  In  both  configurations  particle  contacts 
have  broken,  and  new  contacts  formed. 

Figure  14  and  Figure  15  show  overall  stress-strain  curves  for  this  example.  The 
58-particle  assembly  has  been  subjected  to  a  total  strain  of  cn  =  £22  =  —2.66%  and 
C12  =  £21  =  —3.8%.  Several  observations  already  mentioned  from  the  previous 
examples  may  be  noted  again  from  these  results:  (a)  for  a  constant  H’ Ihx-,  the 
overall  stresses  normalized  with  respect  to  the  same;  (b)  fixing  a  particle  in 

the  assembly  does  not  affect  the  resulting  overall  stress-strain  curve;  (c)  structural 
anisotropy  is  created  by  the  irregular  particle  assembly;  and  (d)  continued  volume¬ 
preserving  shearing  deformation  produces  an  overall  softening  response.  Observation 
(c)  is  evident  in  Figure  14  which  shows  that  during  isotropic  compression,  a\\  does 
not  equal  to  an  even  though  en  equals  £22'  Observation  (d)  may  be  qualified  — 


the  overall  response  may  harden  once  again  as  more  stable  contacts  form  from  a 
collapsed  microstructure,  as  evidenced  by  the  hardening  response  represented  by 
the  tail  of  the  stress-strain  curve  of  Figure  15. 

Finally,  the  overall  stress-strain  curve  of  Figure  16  represents  the  response 
of  the  same  58-particle  assembly  to  the  following  imposed  strain  history:  (a) 
total  isotropic  compression  of  en  =  622  =  —1.0  x  10~^;  (b)  simple  shear 

of  Aei2  =  Ae2i  =  —2.2  x  10“^;  (c)  additional  isotropic  compression  of 

Aeii  =  Ae22  =  —1-0  x  10”^;  and  (d)  simple  shear  of  Aei2  =  Ae2i  =  —2.8  x  10“^. 
The  simple  shear  is  applied  in  steps  of  Aei2  =  Ae2i  =  —1.0  x  10  Both 
the  initial  and  intermediate  isotropic  compression  stages  resulted  in  a  change  in  the 
overall  shear  stress  d’i2  due  to  anisotropy  effects.  The  application  of  an  intermediate 
isotropic  compression  was  necessitated  by  the  deteriorating  numerical  conditioning 
of  the  problem,  manifested  in  the  form  of  lack  of  quadratic  convergence  in  Newton 
iteration,  as  the  overall  stress  response  reaches  a  plateau.  The  additional  isotropic 
compression  is  seen  to  have  resulted  in  a  gain  of  shear  strength. 

One  hundred  ninety  six-particle  random  assembly 

The  assembly  is  composed  of  196  randomly  arranged  circular  disks  of  varying 
sizes  having  a  mean  radius  of  rave  =  0.0092  units  and  contained  in  a  cell  of 
dimensions  0.2566  units  x  0.2564  units.  The  initial  positions  of  the  particles  are 
shown  in  Figure  17,  and  are  identical  to  those  used  in  [36].  The  control  volume  has 
410  initial  contacts  and  an  initial  void  ratio  of  0.1426  (i.e.,  ratio  between  total  area 
of  voids  to  total  area  of  circles)  representing  a  dense  packing.  To  ensure  that  there 
are  no  “numerical  slacks”  between  adjacent  particles  due  to  initial  placement,  and 
that  the  particles  are  indeed  touching  at  the  contact  points,  the  control  volume  was 
isotropically  compressed  to  initial  macroscopic  strains  of  cn  =  622  =  —1.0%. 

The  numerical  algorithm  is  next  tested  for  convergence.  The  control  volume  is 
compressed  further  to  additional  isotropic  strains  of  Aen  =  Ae22  =  —1%  applied 
in  one,  two,  four,  10,  20,  50,  and  100  increments.  Table  2  shows  the  predicted 
incremental  macroscopic  normal  stresses,  Aan  and  Ao'22,  at  cumulative  incremental 
normal  strains  of  Aen  =  ^£22  =  —0.5%  and  Aen  =  ^^22  =  —1.0%.  Note 
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Table  2.  Convergence  test  for  the  196-particle  assembly:  isotropic  compression 
=  kx  =  I  X  10^  units). 


Cumulative  incremental  strain,  Aen  =  Ae22  =  —0.5%: 


Strain  increment,  % 

Stress  Acrii 

Stress  A(T22 

1.00 

n/a 

n/a 

0.50 

45.2789 

40.8908 

0.25 

44.9810 

40.6005 

0.10 

44.7728 

40.4470 

0.05 

44.7013 

40.3923 

0.02 

44.6572 

40.3589 

0.01 

44.6420 

40.3467 

Cumulative  incremental  strain,  Aen  =  Ae22  = 

-1.0%: 

Strain  increment,  % 

Stress  A(t\\ 

Stress  Ad-22 

1.00 

110.8596 

100.1927 

0.50 

109.6300 

99.2609 

0.25 

108.9399 

98.7138 

0.10 

108.4900 

98.4174 

0.05 

108.3350 

98.3125 

0.02 

108.2407 

98.2502 

0.01 

108.2088 

98.2289 

n/a  =  not  applicable 


in  Table  2  that  the  resulting  normal  stresses  are  not  the  same  due  to  anisotropy 
effects.  In  fact,  a  non-zero  macroscopic  shear  stress  of  Aai2  is  also  produced  by  this 
simple  isotropic  compression.  The  results  shown  in  Table  2  demonstrates  that  the 
algorithm  is  convergent  under  an  isotropic  strain  field,  and  that  for  a  volumetric 
compression  of  Aly  —  Aen  +  Ae22  =  —2.0%  the  error  of  the  one-step  solution 
for  the  normal  stress  sum,  Aan  -f  Ao'22,  is  in  the  order  of  2.2%  relative  to  the 
100-step  solution. 

Next,  the  algorithm  is  tested  for  convergence  under  a  simple  shear  strain  field. 
The  initial  condition  for  this  experiment  is  the  final  configuration  that  resulted 
from  the  one-step  solution  of  Table  2  (i.e.,  under  cumulative  normal  strains  of 
cii  =  ^22  =  —2.0%).  From  this  initial  condition,  the  control  volume  is  then  subjected 


Table  3.  Convergence  test  for  the  196-particle  assembly:  simple  shear 
{kjf  =  kx  =  1  X  10'^  units). 


Cumulative  incremental  strain,  Aei2  =  Ae2i  =  —0.5%: 


Strain  increment,  % 

Stress  A<ti2  =  A<T2i 

-0.50 

-21.1606 

-0.25 

-21.1393 

-0.20 

n/a 

-0.10 

-21.1180 

-0.05 

-21.1104 

-0.01 

-21.1028 

Cumulative  incremental  strain,  Aei2  = 

Ae2i  =  -1.0%: 

Strain  increment,  % 

Stress  Aa'12  =  Aa2i 

-0.50 

n/c 

-0.25 

n/c 

-0.20 

-40.0290 

-0.10 

-40.0154 

-0.05 

-40.0078 

-0.01 

-40.0002 

n/a  =  not  applicable;  n/c  =  no  convergence 

to  a  total  simple  shear  strain  of  ei2  =  —1.0%  applied  in  two,  four,  five,  10,  20,  and 
100  steps.  Results  of  the  convergence  study  are  shown  in  Table  3  and  suggest  that, 
again,  the  algorithm  is  convergent  in  the  sense  that  there  exists  a  shear  stress  to 
which  the  solution  tends  as  the  number  of  steps  is  increased.  However,  iterations 
for  the  two-  and  the  four-step  solutions  failed  to  provide  a  convergent  solution  on 
the  last  load  step. 

Note  that  the  strains  in  the  order  of  1%  are  generally  considered  “small”  in  the 
macroscopic  sense.  However,  in  the  microscopic  level,  imposed  macroscopic  strains 
of  this  order  could  produce  significant  relative  motion  between  adjacent  particles, 
and  thus,  can  hardly  be  considered  “small.”  Existing  numerical  solutions  [19-22, 
26-28]  often  limit  the  macroscopic  strain  increments  to  values  well  below  1%  in  order 
that  the  difficulty  associated  with  numerical  conditioning  may  be  circumvented.  An 
important  feature  of  the  proposed  model  lies  in  its  ability  to  provide  convergent 


results  even  for  macroscopic  strain  fields  in  the  order  of  0.1  to  1.0%. 

Using  the  same  initial  configuration  of  Figure  17,  the  particle  assembly  now 
experiences  a  complex  strain  path  characterized  by  alternating  isotropic  compression 
and  simple  shear.  Figure  18  shows  a  plot  of  the  overall  shear  stress  ai2  versus 
overall  tensorial  shear  strain  ei2  for  the  196-particle  assembly.  The  strain  history 
for  this  problem  is  such  that  the  assembly  is  compressed  isotropically  by  an  amount 
of  Aeii  =  Ae22  =  —1.0%  each  time  that  the  shear  stress-shear  strain  curve 

reaches  a  plateau.  Thus,  re-hardening  takes  place  immediately  after  each  isotropic 
compression.  The  deformed  configuration  at  total  strains  of  en  =  £22  =  —3%  and 
£12  =  —4.5%  are  shown  in  Figure  19. 

In  a  separate  parametric  study.  Figure  20  shows  that  the  macroscopic  shear 
stress-shear  strain  response  is  not  significantly  influenced  by  the  assumption  of 
when  a  contact  is  initially  formed  between  two  previously  non-overlapping  particles. 
Recall  that  for  lOPT  =  1,  initial  contact  is  assumed  to  take  place  at  time  in+i  if 
contact  is  detected  between  two  previously  separated  particles;  if  lOPT  =  0,  contact 
is  assumed  to  take  place  at  time  (see  Box  1).  Since  Figure  20  shows  that  the 
predicted  results  are  nearly  the  same,  either  assumption  may  be  employed  in  the 
analysis.  The  use  of  lOPT  =  0  seems  to  generally  lead  to  a  more  stable  numerical 
solution. 


Figure  7.  Initial  configuration  for  the  64-particle 
control  volume;  particle  radius  =  1.0  units. 
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Figure  8.  Isolated  clusters  of  particles  form  from 
an  initially  loose  packing  of  circular  disks;  particle 
radius  =  1.0  units. 
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Figure  9.  Isotropic  compression  followed  by 
lateral  shearing  for  the  64-particle  assembly; 
particle  radius  =  1.0  units. 
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NORMALIZED  SHEAR  STRESS  x  E-4 


SHEAR  STRAIN,  % 


Figure  10.  Normalized  overall  shear  stress 
versus  overall  shear  strain  eu  showing  the  effects 
of  the  moduli  ratio  =  H'  {kx- 
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Figure  11.  Periodicity  test:  deformed  configura¬ 
tion  for  the  60-particle  control  volume;  particle 
radius  =  1.0  units. 
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Figure  12.  Initial  configuration  for  the  58-particle 
control  volume;  particle  radius  =  1.0  units. 
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Figure  13.  Deformed  configuration  for  the 
58-particle  control  volume  showing  zero-energy 
modes  when  the  assembly  is  fully  unconstrained; 
particle  radius  =  1.0  units. 
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NORMALIZED  NORMAL  STRESS  X  E-4 


DIRECTION  11 


0.5  1  1.5  2 

NORMAL  STRAIN,  % 


2.5 


Figure  14.  Normalized  overall  normal  stresses 
aJi/fcj’  and  0-22 /A:r  versus  overall  normal  strains 
eii  and  622,  respectively,  showing  the  effects  of 
structural  anisotropy  for  the  58-particle  assembly. 
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Figure  15.  Normalized  overall  shear  stress 
versus  overall  shear  strain  ei2  for  the  58-particle 
assembly. 
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Figure  16.  Normalized  overall  shear  stress  al^IkT 
versus  overall  shear  strain  ei2  for  the  58-particle 
assembly. 
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Figure  17.  Initial  configuration  for  the  196- 
particle  assembly.  Coefficient  of  uniformity  Cy,  = 
1.71;  coefficient  of  curvature  Cc  =  0.90;  mean 
particle  radius  =  0.0092  units. 


56 


Figure  18.  Overall  shear  stress  d’l2  versus  overall 
shear  strain  ei2  for  the  196-particle  assembly 
subjected  to  alternating  isotropic  compression 
and  simple  shear  strain  history. 


Figure  19.  Deformed  configuration  for  the 
196-particle  assembly;  mean  particle  radius  = 
0.0092  units. 


Figure  20.  Overall  shear  stress  0-^2  versus  overall 
shear  strain  ei2  for  the  196-particle  assembly; 
IOPT=0,  contact  is  assumed  formed  ai  t  =  in 
between  two  previously  non-overlapping  particles; 
IOPT=l,  contact  is  formed  at  t  =  tn+i- 


Chapter  5 

Stress-driven  Problem 


5.1.  Introduction 

When  a  macroscopic  point  in  a  granular  material  experiences  an  overall  stress 
history,  it  responds  according  with  some  uniform  macroscopic  deformation.  This 
overall  motion  must  be  compatible  with  the  movement  of  the  individual  particles 
contained  in  the  control  volume.  Thus,  when  given  a  stress  history,  the  macroscopic 
uniform  deformation  of  the  particle  assembly  and  the  discrete  microscopic  particle 
motions  constitute  the  principa.1  unknowns.  In  the  language  of  computational 
plasticity,  one  has  a  stress- driven,  or  inverse,  problem.  A  stress-driven  problem 
requires  two  levels  of  analysis: 

i.  micromechanical  level  to  calculate  the  contact  forces  and,  eventually,  the 
overall  stresses  &  resulting  from  an  imposed  overall  deformation  gradient  A, 

ii.  macromechanical  level  to  iteratively  determine  the  uniform  control  volume 
motion  A*  which  exactly  produces  the  prescribed  overall  stresses  a*. 

The  strain-driven  algorithm  presented  in  Chapter  4  forms  the  basis  of  the 
micromechanical  level  by  determining  the  particle  contact  forces  for  a  given  overall 
deformation  of  the  control  volume. 

It  is  emphasized  that  while  the  formulation  in  this  chapter  admits  finite  motions 
of  the  individual  granules,  the  overall  material  response  is  based  on  infinitesimal 
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theory,  i.e.,  A  is  assumed  to  be  ‘small.’  Consequently,  the  present  formulation 
captures  only  the  material  overall  response.  Finite  rotation  and/or  large  deformation 
of  the  control  volume,  which  could  be  important  in  many  boundary-value  problems, 
may  be  readily  incorporated  in  the  macroscopic  formulation  in  a  future  extension 
of  the  model. 

This  chapter  presents  alternative  algorithms  for  solving  the  inverse  problem 
based  on  Newton  and  Newton- type  methods.  First,  a  description  of  the  algorithm  for 
the  solution  of  the  stress-driven  problem  is  presented.  The  remainder  of  the  chapter 
considers  the  issues  surrounding  the  consistent  linearization  of  the  overall  elasto- 
plastic  constitutive  relations.  In  order  to  exactly  derive  the  algorithmic  tangent 
operator,  an  additive  decomposition  of  the  particle  motions  is  presented.  The 
implications  of  this  decomposition  on  particle  contact  slip,  contact  forces,  and  the 
overall  stress  are  then  investigated  to  render  the  derivation  of  the  exact  algorithmic 
tangent.  This  chapter  then  devotes  special  attention  to  the  implementational  hurdles 
surrounding  the  exact  tangent  operator.  Finally,  to  circumvent  these  difficulties,  a 
secant  approximation  is  described  that  retains  the  essential  properties  of  the  exact 
description. 


5.2.  Algorithm  for  problems  with  imposed  overall  stresses 

Let  be  the  imposed  overall  stresses.  For  equilibrium,  one  can  write  the 

macroscopic  stress  equation  at  time  station  as 

=0.  (5-1) 

where  denotes  the  overall  displacement  gradients  compatible  with  the  imposed 

overall  stresses  and  V’n-l-i  represents  a  set  of  strictly  micromechanical 

variables.  This  equilibrium  equation  represents  an  extension  of  the  functional 
relation  (3.29)  to  include  the  effects  of  the  a.ssembly’s  microstructure.  The  physical 
meaning  of  equation  (5.1)  is  that  the  applied  loads  equilibrate  the  internal  stress 
which  depends  in  a  very  complex  manner  upon  the  uniform  control  volume  motion 
and  upon  micromechanical  variables  such  as  the  particle  coordination  numbers,  the 


formation  and  rupture  of  particle  contacts,  and  the  stability  of  the  control  volume’s 
microstructure. 

In  plane  strain  applications,  the  applied  stress  represents  a  vector  with  four 
independent  components.  However,  since  the  balance  of  angular  momentum  (3.4) 
requires  that  the  calculated  overall  stress  be  symmetric,  equation  (5.1)  yields  three 
coupled  nonlinear  equations  for  each  component  of  the  applied  stress: 


^n+1  —  {<^115  <^225  ^12}n+l  •  (^-2) 

The  unknown  vector  of  the  uniform  motion  then  becomes 

^n+l  =  {All,  A22,  Ai2}^+1  .  (5.3) 

The  applied  stress  vector  in  (5.2)  will  be  formed  by  instantaneously  imposing  the 
incremental  stress  ^j^^i  —  at  discrete  time  steps  At  =  tn+i  —  tn-  Likewise,  the 
internal  stress  &  instantaneously  varies  following  the  application  of  the  incremental 
stresses. 

We  now  denote  {GEXT)n+l  as  the  vector  of  applied  overall  stresses  given  by 

{GEXT)n+l  -  ^n+l  ,  (5.4) 

and  GjffT  as  the  vector  of  internal  stresses  given  by 

Gint  =  ^(A„+i,  V’n+i)  •  (5.5) 

The  condition  of  macroscopic  stress  equilibrium  can  then  be  rewritten  in 
residual  form  in  terms  of  macroscopic  quantities  as 

R  =  Gj^riAl+i)  -  {GEXT)n+i  =  0  (5.6) 

where  are  the  roots  of  the  nonlinear  system  of  equations  (5.1).  Linearizing 

the  residual  equation  gives 


=  {GEXT)n+l  -  Gixt{AI^i) 


(5.7) 


where 


gW(4+i)  =  c  = 

is  the  tangent  matrix,  and  AA^  represents  the  A:th  search  direction  for  the  solution 
vector  Note  that  the  matrix  C  in  (5.8)  represents  the  material  stress-strain 

matrix  when  the  elements  of  are  small. 


Applying  Newton’s  method,  the  next  estimate  of  can  be  computed  from 

=  ^‘+1  -  -K*-  (S-9) 


The  method  is  said  to  have  converged  when 


||fi‘||/||il"||<RTOL. 


(6.10) 


where  RTOL  is  a  prescribed  error  tolerance. 

Box  3  summarizes  the  necessary  steps  for  solving  a  stress-driven  problem.  In 
view  of  the  two  levels  of  analysis  required  to  solve  this  type  of  problem,  two  nested 
Newton  iteration  loops  are  necessary.  Inclusion  of  the  call  to  the  micromechanical 
level  clarifies  the  central  role  played  by  the  strain-driven  algorithm.  Note  in  Box  3 
that  the  inverse  algorithm  is  only  at  best  as  stable  as  the  inner  Newton  loop  for  the 
strain- driven  algorithm. 


5.3.  Particle  displacement  decomposition 

Consider  the  exact  evaluation  of  the  algorithmic  tangent  (5.8).  Central  to 
this  computation  lies  the  complex  micromechanical-macromechanical  connection 
between  the  independent  particle  motions  and  the  overall  control  volume 
deformation.  The  micromechanical  level  subjects  the  assembly  of  particles  to 
a  prescribed  macroscopic  uniform  motion  and  then  allows  the  particles  to  move 
in  a  microscopic  sense  to  achieve  equilibrium  [26-28].  This  section  derives  a 
decomposition  of  the  particle  displacements  to  formalize  this  micromechanical- 
macromechanical  connection  essential  for  exactly  evaluating  the  algorithmic  tangent. 


Box  3.  Stress- driven  algorithm. 


Macromechanical  level 

1.  Given  or*,  and  k  =  0. 

2.  Compute:  G'j . 

^n+i  =  ^n+1  ~  [^/jvrl-^n-l-l)]  r^^- 

4.  Call  Micromechanical  level  in  Box  2  with  An,  and  i  =  0. 

5-  ^  =  5E"=.(r®i“+i"®n- 

6.  il‘+>  =  o-(Aj+.)  - 

7.  IF  ||i2^+i||/||J7°||  <  RTOL,  RETURN. 

8.  k  ^k  +  1  and  GO  TO  2.  □ 


To  this  end,  the  relationship  between  the  particle  displacements  and  the  overall 
displacement  gradient  must  first  be  constructed. 

Recall  that  the  dependence  of  the  contact  forces  on  the  imposed 
displacement  gradient  A  followed  directly  from  the  initial  imposition  of  the  uniform 
motion  on  the  particle  centroids.  This  motion  causes  contacting  particles  to  either 
slip,  overlap,  or  separate  creating  unbalanced  contact  forces  and  acts  as  the  initial 
estimate  for  the  Newton  solution  given  in  Box  2.  The  perturbed  momentum  balance 
equations  given  in  (4.13)  as 

r  =  FiNT{A-X^)^0  (5.11) 

create  the  residual  in  the  force  vector  which  must  be  dissipated  by  the  Newton 
iteration.  Thus,  the  initial  particle  displacements  (i.e.,  the  first  guess  of  the  new 
particle  equilibrium  positions)  depend  solely  on  the  uniform  motion  A. 

The  continued  dependence  of  on  the  displacement  gradient  A  follows  from 
the  dissipation  of  r  in  (5.11)  through  the  Newton  iteration  as  the  particles  displace  to 
the  new  equilibrium  configuration.  The  sides  of  the  unit  cell  remain  fixed  according 


to  the  uniform  motion  of  A.  In  other  words,  after  the  initial  imposition  of  A,  the 
particle  displacements  to  equilibrium  continue  to  depend  on  A,  but  they  also  depend 
on  the  local  effects  of  the  material’s  microstructure  described  by  the  variables 
Thus,  from  a  typical  particle’s  perspective,  the  dependency  on  the  displacement 
gradient  A  and  the  micromechanical  variables  can  be  formalized  as 

dn-f-l  —  .  (5.12) 

This  displacement  relationship  can  be  understood  from  a  different  perspective. 
Consider  a  densely  packed  control  volume  experiencing  a  small  overall  shearing 
motion  that  slightly  perturbs  the  momentum  balance  equations  (5.11).  The 
constituent  particles  need  only  displace  and  rotate  a  small  amount  depending  on 
the  current  set  of  micromechanical  variables  if)  to  reach  equilibrium.  Therefore,  in 
the  limit  as  the  magnitude  of  the  overall  motions  approaches  a  small  value,  say  e,  the 
uniform  particle  motions  may  accurately  describe  the  total  particle  displacements, 

lim  dn-i-i  =  dn+i{An+i)  ■  (5.13) 

IMn+lIHe 

Since  the  macroscopic  motion  is  not  restricted  to  be  vanishingly  small,  the  particle 
displacements  will  also  depend  on  the  microscopic  variables  if),  and  therefore 

dji+l  —  dn-^l{An-{-l,  lf)n+lf  • 

This  displacement  relationship  can  be  made  explicit  by  additively  decomposing 
the  particles’  displacement  degrees  of  freedom  into  a  displacement  d  due  to  the 
uniform  control  volume  motion  and  a  microscopic  displacement  d  describing  the 
particles’  equilibrium  motion 

d^{An+l,‘tkn+l)  =  d^^An+l)  +  d^{An-\-\^  '0n+l)  5  (5.15) 

where  particle  A’s  point  of  view  has  been  taken.  Note  that  the  microscopic 
displacement  d  also  depends  on  the  uniform  motion. 

In  two- dimensions,  the  associated  displacements  for  particles  A  and  B 
connected  by  contact  element  ‘e’  can  be  written  as 

(f  =  d^  +  d^,  (5.16) 
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where 


^(6x1) 


~d^  1 

d^  1 

d^ 

^  ’  ^^(6x1)  -  ■* 

0 

~dB 

A 

^  1  ^(6x1)  ~  ^ 

d^ 

.  ^n+1  . 

.  0  > 

^  ^n+1  > 

(5.17) 


Note  that  according  to  this  decomposition  the  uniform  control  volume  motion  results 
in  only  particle  translation.  The  individual  particle  rotations  result  from  the  locally 
discrete  particle  motion  to  reach  equilibrium. 


5.4.  Displacement  decomposition  and  particle  kinematics/contact  forces 

This  section  considers  the  implications  of  the  particle  displacement 
decomposition  for  the  particle  contact  slip  and  the  contact  forces.  It  pays  special 
attention  to  the  impact  of  the  decomposition  on  the  path  of  the  problem  at  a  typical 
particle  contact. 

Consider  application  of  the  displacement  decomposition  (5.15)  to  the  definition 
of  the  contact  unit  normal.  As  seen  in  Chapter  2,  the  definition  of  the  contact  normal 
is  central  to  the  kinematic  description  of  particle  contact  slip.  Again  taking  particle 
A’s  point  of  view,  let  the  contact  normal  in  two  dimensions  be 

^n+\  =  =  Wl/|l^n+l||  =  (5.18) 

where 

=  +  +  (5.19) 

=  X®  +  -  (X.*  +  ) .  (5.20) 

Similarly,  the  contact  normal  after  the  application  of  the  uniform  macroscopic 
motion  is  given  by 


(5.21) 


where 


h  =  x^-)-dS+i-{x''  +  di+,). 


(5.22) 


The  normal  indentation  6  can  be  written  as 

(5„+i=r"^  +  r^-||Z„+l||,  (5.23) 

where  and  denote  the  radius  of  particle  A  and  B  respectively. 

The  definitions  of  and  nd  may  be  used  to  segregate  the  total  slip 

increment  A7  into  slip  due  to  the  uniform  control  volume  motion  Ajm  and  slip 
due  to  the  particle  equilibrium  motion  A'yf^ 


A7  =  A^n  + 

(5.24) 

Recall  that  the  total  slip  increment  is  written  as 

A7  =  A0^^r^  +  AO^^r^  , 

(5.25) 

where 

(5.26) 

A»»^  =  Ae'^-(Ci-0> 

(5.27) 

and 

Ae^  =  sin  ^(63  •  x  n;^^.i) . 

(5.28) 

Incorporation  of  and  nd  into  these  equations  leads  naturally  to  the  following 
definitions: 


A^n  =  Aefr^  +  Ae?^rB, 
AOf  =  AC  , 

Aei^  =  Ae^ , 

AO^  =  sin“^(e3  •  n„  x  na) , 


A7*  =  A^d^r^  +  Ae?^r^  ; 

Aef  =  A^f  -  {0^^,  -  ei) ; 
A0f^  =  A0?-(Ci-^f); 

A6?  =  sin“^(e3  •  x  n„+i) 


(5.29) 


It  is  important  to  note  that  equations  (5.24)  and  (5.29)  only  redefine  the 
definition  of  the  total  slip  A7  in  terms  of  the  displacement  decomposition.  The  path 
of  the  problem  remains  unaltered  in  the  sense  that  the  solution  moves  from  time 


station  tn  to  tn+i  without  updating  the  path  dependent  variables  at  an  intermediate 
station,  say,  For  example,  the  contact  plastic  slip  is  defined  for  the  entire 
increment  At  =  tn+i  —  tn  as 


1A7^1  = 


\f^\  -a„  +  /jvtan?i 
Hekr  +  H' 


where 

f^  =  {fT)n  +  H,kTAr 


(5.30) 


(5.31) 


The  loading/unloading  conditions  have  not  been  used  to  compute  intermediate 
values  such  as  A7?  and  A7?.  Using  the  displacement  decomposition  as  the  basis 
for  intermediate  nonconverged  states  appears  questionable  for  a  problem  which, 
physically,  is  path-dependent  [37].  If  unloading  of  a  contact  element  occurs  during 
the  Newton  iteration  process,  each  new  iteration  should  start  from  the  converged 
configuration  dn  and  not  from  an  intermediate  configuration  like,  say,  da- 

Likewise,  it  is  important  to  note  that  loading/unloading  conditions  are  not 
employed  to  derive  an  intermediate  force  resulting  from  the  uniform  motion,  say 
(/j’)ii,  and  a  remainder  force,  say  {fT)h-  The  displacement  decomposition  does  not 
result  in  a  decomposition  in  the  contact  forces  such  as 


{fT)n+l  {fT)h{d)  -f  (/r)n(d,  d).  (5.32) 


Such  an  application  of  the  displacement  decomposition  would  erroneously  change 
the  path  of  the  problem  from 


tn 


to  a  new  path  described  by 


(5.33) 


(5.34) 


where 


means  that  the  force  fn+i  is  updated  from  configuration  at  tn  to  tn+i- 


5.5.  Displacement  decomposition  and  the  macroscopic  stresses 

Now  consider  the  implications  of  the  displacement  decomposition  on  the 
calculated  stress.  In  Section  3.4,  the  initial  description  of  the  internal  stress  was 
simply 

^n+i  =  a-n-k-i{A).  (5.35) 

With  the  inclusion  of  the  micromechanical  variables  tj),  this  expression  then  became 

d-n+i  =  d-„+i(A,  V’)  •  (5.36) 

A  thorough  examination  of  the  overall  stress  from  a  different  perspective  leads  to 
a  more  precise  expression  for  the  dependency  of  the  stress.  The  expression  for  the 
stress  (3.12)  shows  a  dependency  on  the  particle  contact  forces  and  the  location  of 
the  particle  centroid, 

a  =  (5.37) 

With  the  displacement  decomposition,  the  contact  forces  T  may  be  written  as 

T^T(d,d).  (5.38) 

Thus,  the  calculated  stress  can  be  then  viewed  as  a  fundamental  function  of  the 
uniform  and  equilibrium  displacements 

&n+i  =  ^n+i{d,  d) .  (5.39) 

This  functional  relationship  can  be  reconciled  with  equation  (5.36)  by  substituting 
the  displacement  decomposition  to  yield 

^n+i  =  d-n+i{d{A),  d{A,  V*)) .  (5.40) 

Note  that  equation  (5.38)  can  be  easily  demonstrated  by  considering  the 
expression  for  the  tangential  contact  force 


(/r)n+i  =  (/r).  +  H.kr  (A7  -  A^’’) . 


(5.41) 


Substituting  (5.24),  (5.30),  and  (5.31)  yields 
Hekr 


(/r)n+i  = 


H,kT  +  H'  L 
+ 


H'iA'jn  -  A7ji)  -  sign(-)a„  +  sign{-)f^  tan  (j) 
H' 


(5.42) 


(/r)n 


H,kT  +  H' 

where  (•)  =  (/r)n  +  H^kj'A'f.  From  (5.29),  it  follows  that  A7ji  =  Ajn{d)  and 
A7n  —  A7ji(d,  d),  and  therefore 

{fT)n+i  =  fTid,d)  (5.43) 

which  leads  directly  to  equation  (5.38). 


5.6.  Exact  algorithmic  tangent  operator  and  concentration  tensor 

Now,  with  the  expression  (5.39)  for  the  stress,  an  exact  expression  for  the 
tangent  operator  can  be  derived.  Employing  the  chain  rule  on  the  calculated  stress 
&  gives  (dropping  the  subscripts  n  +  1  and  superscripts  k) 


,  f  A  \  _  do-{d,d)  da  dd  da  dd 

Gwr(A)  =  +  iSal  ■ 

The  derivatives  of  the  stress  a  can  be  explicitly  derived  from  (3.12)  as 


(5.44) 


ddk 

da, 


ii  =  f 4.  J.  _L 

9dk  '  ^  \ddk  '  '  ddk  ddk  ‘  ddk 


(5.46) 


Appendix  D  contains  a  complete  discussion  of  these  force  derivatives. 

The  displacement  derivatives  in  equation  (5.44)  can  be  summarized  in  terms 
of  the  displacement  decomposition  as 

dA  dA  dA 


(5.47) 


Taking  the  point  of  view  of  particle  A,  the  first  derivative  on  the  right-hand  side 
can  be  found  from 


d^  =  A-Xt 


(5.48) 


The  derivative  can  be  written  as 

ddA 

dA 


'Xf  0 
0 

0  0  0 


(5.49) 


Equation  (5.47)  represents  the  change  in  the  total  particle  motion  with 
respect  to  the  change  in  the  overall  motion  of  the  control  volume  —  the  so-called 
‘concentration  tensor’  [23-28,38]. 

The  exact  evaluation  of  the  concentration  tensor  requires  explicitly  extending 
the  displacement  decomposition  into  the  micromechanical  level  residual  r  which  can 
be  written  as 


r  =  F/jvT(<^n-t-i)  =  0  (5.50) 

=  FjNT{d,d)  (5.51) 

^el 

=  \Jf{d/d).  (5.52) 

e=l 

The  first  variation  of  this  expression  yields 

^8d-\-  ^r8d  =  0.  (5.53) 

dd  dd 


where  8  denotes  a  small  variation.  Since  the  displacements  d  and  d  depend  on  the 
macroscopic  motion  A,  use  of  the  chain  rule  on  (5.53)  gives 


dr  dd  dr  dd 

dd  9 A  dd  9 A 


(5.54) 


It  follows  that  for  an  arbitrary  8  A 


dr  dd  dr  dd 
dd'^  dd  dA 


(5.55) 


Substituting  equation  (5.52)  into  (5.55)  yields  an  exact  expression  for  the 
concentration  tensor  as 

^  =  (5.56) 

where 

71^1 


^=U*'=U 


e—1 

riel 


dr 

^1 


K  +  K 


=  U(*‘  +  *')  =  U 

e=l  e=l 


dr 

dd^ 


(5.57) 

(5.58) 


The  element  contributions  to  the  matrices  K  and  K  can  be  found  in  Appendix  D. 
The  definition  of  the  algorithmic  tangent  matrix  which  can  now  be  written  as 


da  da 
^~~dd 


f)d\  J 


dd 

dA 


(5.59) 


As  seen  from  Box  3,  the  algorithmic  tangent  will  be  evaluated  after  the 
micromechanical  level  has  reached  convergence  for  the  given  macroscopic  motion. 

The  numerical  simulations  of  Section  4.4  showed  that  when  a  control  volume 
with  an  initially  collapsible  structure  experiences  a  shearing  deformation  beyond 
its  ability  to  compact,  particles  separated  into  isolated  clusters.  Numerically,  this 
phenomenon  manifested  itself  as  zeros  on  the  diagonal  of  the  factorized  microscopic 
tangent  operator.  These  zeros  are  not  operated  upon  during  the  back  substitution 
process.  The  stress-driven  problem  does  not  preclude  such  a  phenomenon.  In 
fact,  such  a  ‘strain  softening’  effect  could  have  disastrous  implications  for  the  exact 
evaluation  of  the  tangent  and  in  particular  the  evaluation  of  K~^ .  If  K  becomes 
singular,  the  zeros  on  the  diagonal  cannot  be  ignored  and  will  make  the  calculation 
of  the  exact  algorithmic  tangent  intractable. 

In  addition,  the  complexity  of  the  exact  algorithmic  tangent  should  be  noted. 
Equation  (5.59)  requires  the  considerable  effort  to  evaluate  the  derivatives  contained 
in  K  and  K  and  also  the  explicit  computation  of  K~^.  The  matrices  fe®  and  fe® 
represent  contact  element  contributions  and  therefore  must  be  computed  for  every 
particle  contact  in  the  control  volume.  Therefore,  the  intensive  computational  and 


assembly  effort  expended  in  the  exact  evaluation  of  K,  and  K  also  represents 

a  significant  difficulty  engendered  with  the  implementation  of  the  exact  algorithmic 
tangent  operator. 


5.7.  Discretized  Newton  method 


Although  Newton’s  method  is  theoretically  attractive  for  the  stress-driven 
problem,  it  is  expensive  and  difficult  (if  not  impossible)  to  implement.  Each  iteration 
not  only  requires  the  calculation  of  the  components  of  the  macroscopic  residual  but 
also  the  components  of  the  algorithmic  tangent.  As  Section  3.5  and  Appendix  D 
show,  the  partial  derivatives  in  the  tangent  matrix  do  not  have  a  simple  functional 
form.  The  most  direct  approach  to  circumventing  the  explicit  computation  of  these 
derivatives  is  simply  to  approximate  the  algorithmic  tangent  by  difference  quotients. 
In  what  follows,  two  such  approximations  for  the  general  stress-driven  problem  in 
three  dimensions  are  considered. 


Two  commonly  used  difference  approximations  to  the  algorithmic  tangent 
operator  take  the  form 

dRiiA)  .  1 


and 


^R^{A)  ^  J_ 

^.A  j  ^ij 


=  -j—  Ri  -f  hije^^  —  Ri  (A)  , 


jfc=i 


j 

Ri  ^  —  Ri  (^A  +  ^  hik^  ^ 

k—1 


(5.60) 


(5.61) 


where  hij  denotes  discretization  parameters,  and  represents  the  jth  coordinate 
vector.  The  residual  R  and  parameter  vector  h  have  the  mapping  C  3?®  and 
h.  C  X  3?®  respectively. 

Let  ^{A,h)  denote  the  difference  approximations  (5.60)  and  (5.61)  with  the 
property  that  whenever  algorithmic  tangent  exists,  then 


lim  ^(A,  h)  =  Grjj^j'i^A,  .  (5.62) 

h—*0 
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(5.63) 


Thus,  the  next  estimate  of  can  be  computed  from 


1 


This  iteration  scheme  is  the  so-called  ‘discretized-Newton’  method  and  represents  a 
special  case  of  the  general  secant  method.  Note  that  the  discretization  parameters 
h  vary  with  each  iteration. 

Now  consider  several  possible  ways  of  choosing  h.  If  hij  =  hi  such  that 

-  A'^  ,  (5.64) 


the  difference  approximation  (5.60)  yields 
dR{A) 


dA^ 


1 


^  [b(.4‘  +  hte^)  -  h(a‘) 


hi 


(5.65) 


and  the  approximation  (5.61)  gives 

dR{A^) 


dA^ 


+  .... 


6  j=i  j=l 


(5.66) 


The  auxiliary  points  A^  +  (A^  ^  —  Aj)e^  and  A^  +  ^  —  A^)e^  in  these 

difference  approximations  denote  the  points  in  3?®  used  to  approximate  the  tangent 

The  methods  defined  by  (5.65)  and  (5.66)  retain  the  essential  properties  of 
Newton’s  method,  and  have  satisfactory  local  convergence  theorems  [34].  An 
important  property  of  the  method  defined  by  (5.63)  is  that  it  exhibits  the  same 
quadratic  convergence  as  Newton’s  method  while  not  requiring  any  derivatives  of 
the  residual  function  R. 

Note  that  for  a  plane  strain  formulation,  the  approximation  defined  by  (5.65) 
requires  three  function  evaluations  of  the  macroscopic  residual  R.  However,  since 


R{A^~^)  is  available  from  the  previous  iteration,  the  approximation  (5.66)  requires 
only  two  new  evaluations  of  R.  Thus,  with  only  two  evaluations  of  the  macroscopic 
residual  i2,  one  can  circumvent  the  intensive  calculations  for  the  exact  algorithmic 
tangent. 


5.8.  Numerical  simulations 

This  section  presents  the  results  of  two-dimensional  plane-strain  simulations 
on  granular  assemblies  composed  of  either  regular  or  random  initial  packing  of 
circular  disks.  It  also  demonstrates  some  fundamental  properties  of  our  numerical 
model  such  as  quadratic  convergence  of  the  discretized-Newton  approximation,  the 
invariance  of  the  model  under  rigid  body  rotation  of  the  control  volume,  and  the 
prediction  of  the  onset  of  localization,  as  well  as  reports  the  numerical  difficulties 
encountered  during  the  simulation  process. 

Unless  otherwise  stated.  Table  4  defines  the  model  parameters  used  in  all 
analyses  for  the  macroscopic  level  and  the  microscopic  level.  The  overall  stresses 
have  been  calculated  using  a  volume  equal  to  the  initial  cross  sectional  area  Aq 
times  a  unit  thickness.  The  error  tolerances  RTOL  and  rtol  used  in  (4.8)  and  (5.10) 
respectively  ensure  that  the.  iteration  has  sufficiently  converged,  thus  rriinimizing 
the  propagation  of  numerical  errors.  For  the  examples  considered  throughout  this 
section,  the  discretized  Newton  approximation  (5.66)  has  been  employed.  The  same 
results  can  be  obtained  with  the  approximation  (5.65). 

Sixty  four-particle  closest  packing  assembly 

The  initial  granular  assembly  configuration  is  shown  in  Figure  21.  In  this 
example,  the  assembly  consists  of  sixty  four  uniform  circular  particles  of  radius 
r  =  1.0  units  arranged  in  a  closest  packed  configuration.  The  same  configuration 
can  be  obtained  by  rotating  the  control  volume  by  60°.  The  particle  on  the  lower  left- 
hand  corner  has  been  fixed  against  translation  and  rotation  to  arrest*  the  zero  energy 
modes  present  in  the  assembly.  The  two  initial  auxiliary  points  A®  and  are 


Table  4.  Model  parameters. 


Macromechanical  level 

Error  tolerance:  RTOL  =  1.0  x  10~^. 
Micromechanical  level 

Error  tolerance:  rtol  =  1.0  x  10~®. 

Normal  spring  stiffness:  kff  =  1.0  x  10'^. 
Tangential  spring  stiffness:  kx  =  1.0  x  10^. 
Particle  friction  angle:  (f>  =  30°. 

Particle  contact  cohesion:  oq  =  0. 

Contact  hardening  parameter:  H'  =  100. 
Ramp  function  parameter:  c  =  0.10. 


arbitrarily  chosen  as  {-0.02,-0.02,0.0}^  and  {—0.0275,-0.0275,-5.5  x  10“®}^, 
respectively,  to  launch  the  discretized  Newton  approximation  (5.66).  The  converged 
algorithmic  tangent  from  the  previous  time  step  is  employed  as  the  initial  estimate 
of  the  tangent  operator  for  each  new  time  step.  Otherwise,  these  examples  follow 
the  Newton  algorithm  outlined  in  Box  3,  i.e.,  update  the  algorithmic  tangent  each 
iteration. 

The  particle  assembly  experiences  a  stress  history  consisting  of  one  increment 
of  isotropic  compression  Ad-Jj  =  Ad22  =  —200  *  Aq  =  —0.9021 ,  Ad^j  = 
Ad’{2  =  0.0 ,  followed  by  35  increments  of  shear  stress  Ad{j  = 
Ad22  =  0.0;  Ad^i  =  Ad22  =  — 10.0*  Aq  = —0.0451.  The  plot  of  the  normalized 
overall  shear  stress  d{2/^r  versus  the  uniform  shear  motion  A12  appears  in  Figure  22 
and  shows  the  stress  point  reaching  a  yield  plateau.  At  this  plateau,  the  localization 
function  (3.41)  becomes  negative  at  the  last  time  step  indicating  that  the  localization 
criterion  is  satisfied  at  an  intermediate  point  between  the  last  two  time  steps.  The 


deformed  configuration  at  the  end  of  the  loading  history  appears  in  Figure  23. 

The  same  analysis  was  rerun  with  a  shear  stress  loading  in  the  opposite 
direction:  Aal2  =  10.0  *i4o  =  0.0451.  The  same  results  for  the  final  configuration, 
contact  forces,  and  uniform  motion  were  obtained  under  the  appropriate  rotation. 
Thus,  the  model  is  invariant  under  rigid  body  rotations  and  can  capture  the  isotropy 
of  the  assembly  packing. 

Fifty  eight-particle  irregular  assembly 

The  initial  control  volume  configuration  for  this  example  is  shown  in  Figure  12. 
This  control  volume  represents  a  first  step  towards  the  general  analysis  of  randomly 
arranged  and  randomly  sized  particles.  This  repeating  cell  configuration  isolates 
the  effect  of  a  random  arrangement  of  particles  on  the  overall  response  of  the 
material  and  the  performance  of  our  model.  All  particles  in  the  control  volume 
have  a  radius  r  =  1.0  units.  For  these  experiments,  the  two  auxiliary  points 
are  arbitrarily  chosen  to  be  and  to  be  {-0.002,-0.002,  —1.0  x  10“^}  and 
{-0.00175,  -0.00174,  -2.55  x  10"®  ,  respectively. 

For  the  examples  considered  here,  the  control  volume  experiences  six 
increments  of  isotropic  compression.  The  first  two  increments  consist  of  = 

Aa22  =  —200.0  *  Aq  =  —0.7842,  and  the  next  four  increments  consist  of 

Aajfj  =  Ad-22  =  —900.0  *  Ao  =  —3.5291.  The  isotropic  compression  produces 
prestressing  effects  on  the  elastic  springs  and  prevents  the  particle  contacts  from 
breaking  during  the  shearing  process.  The  compression  may  be  thought  of  as 
stabilizing  the  particles  which  may  initially  lie  in  an  unstable  configuration. 

The  particle  assembly  now  experiences  a  shear  stress  loading  in  increments  of 
Ad-12  =  —40.0  *  Ao  =  —0.1568  for  different  values  of  the  hardening  parameter 
H' .  Figure  24  shows  the  variation  of  the  normalized  shear  stress  aj‘2/A:r  with  the 
shear  motion  An  for  each  choice  of  H'.  The  straight  line  portions  of  each  response 
represent  the  elastic  stretching  of  the  tangential  springs.  As  H'  becomes  smaller, 
the  overall  response  approaches  an  elastic-perfectly  plastic  behavior. 

With  the  same  initial  configuration  and  a  hardening  parameter  H'  =  100.0,  the 
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control  volume  is  subjected  to  a  cyclic  shear  loading  consisting  of  80  increments  of 
Ad‘j2  =  — 10.0*Ao  =  —0.0392;  followed  by  180  increments  of  Ad-J2  =  10.0  *Ao  = 
0.0392;  and  finally  191  increments  consisting  of  Ad-J2  =  —10.0  *  Ao  =  —0.0392. 
Figure  25  shows  a  plot  of  the  normalized  overall  shear  stress  response  versus  the 
uniform  shear  motion  Ai2-  This  overall  response  clearly  indicates  a  hardening 
response  resulting  from  the  incorporation  of  particle  contact  hardening.  Note  that 
localization  does  not  occur  during  any  stage  of  loading,  unloading,  or  reloading.  The 
plot  of  the  localization  function  near  the  region  of  the  function  minima  and  at  the 
end  of  the  first  loading  appears  in  Figure  26. 

One  hundred  ninety  six-particle  random  assembly 

The  initial  positions  of  the  particles  are  shown  in  Figure  17.  For  the  examples 
in  this  section,  the  two  auxiliary  points  are  chosen  to  be  {-0.02,-0.02,0.0}* 
and  {-0.0275,-0.0275,-5.5  x  10"®}*  ,  respectively.  To  ensure  an  initially  stable 
initial  configuration  and  that  the  particles  indeed  touch  at  contact  points  the 
control  volume  was  isotropically  compressed  to  an  initial  macroscopic  stress  of 
AaJi  =  Aa|2  =  -5.0  *Ao  =  -75.9968. 

The  numerical  algorithm  is  next  tested  for  convergence.  Following  the  Newton 
algorithm  in  Box  3,  the  control  volume  is  compressed  with  an  additional  isotropic 
stress  of  AaJj  =  Ad-22  =  — 5.0*.4o  =  —75.9668  applied  in  one,  two,  four,  10,  20,  and 
50  increments.  Table  5  shows  the  predicted  incremental  macroscopic  control  volume 
motions  AAn  and  AA22,  at  cumulative  incremental  normal  stresses  of  Ad-{j  = 
Aa^2  =  -2-5  *Ao  =  -37.9984  and  Aal^  =  Ad-^2  =  *  ^0  =  -75.9968. 

Note  that  the  normal  uniform  motions  of  the  control  volume  are  not  equal  due  to 
anisotropy  effects,  and  a  non-zero  macroscopic  shear  motion  AA12  =  A/i2i  is  also 
produced  by  the  isotropic  stress  history.  The  results  shown  in  Table  5  demonstrate 
that  the  inverse  algorithm  is  convergent  under  an  isotropic  stress  field  and  that  for 
a  volumetric  stress  of  Ad-„  =  Ad-J^  -f  Ad'22  =  —10.0  *  Aq  =  —151.9936  the  error  of 
the  one-step  solution  for  the  normal  motion  sum,  AAu  -f  AA22,  is  in  the  order  of 
1.0%  relative  to  the  50-step  solution. 

Table  5  also  shows  the  convergence  profile  of  the  one-step  solution.  This  profile 


Table  5.  Convergence  test  for  the  196-particle  assembly:  isotropic  compression 
(fcjv  =  fey  =  1  X  10'^  units). 


Cumulative  incremental  stress,  Acru  =  A(T22  =  - 

-2.5  *Ao  =  -37.9984  : 

Stress  increment 

AAii,  % 

AA22,  % 

-5.00 

n/a 

n/a 

-2.50 

-0.31043 

-0.34583 

-1.00 

n/a 

n/a 

-0.50 

-0.31160 

-0.34810 

-0.25 

-0.31176 

-0.34841 

-0.10 

-0.31186 

-0.34859 

Cumulative  incremental  stress,  Aan  =  A<722  =  - 

-5.0  *Aq  =  -75.9968  : 

Stress  increment 

AA„,  % 

AA22,  % 

-5.00 

-0.57729 

-0.64472 

-2.50 

-0.57976 

-0.64837 

-1.00 

-0.58129 

-0.65089 

-0.50 

-0.58184 

-0.65176 

-0.25 

-0.58213 

-0.65220 

-0.10 

-0.58230 

-0.65246 

n/a  =  not  applicable 

Convergence  profile  for  stress  increment,  Atrn  = 

A(T22  =  —5.0  *  Aq  =  —75.9968: 

Iteration 

Residual 

1 

7.0710678 

2 

1.2514259 

3 

1.6524780  X  10“^ 

4 

4.0669135  X  10-^ 

5 

2.0566571  X  10"^ 

6 

3.0604956  x  10”® 

typifies  the  performance  of  the  discretized  NewtOn  method  used  for  all  numerical 
simulations  presented  thus  far.  The  approximation  (5.66)  to  the  algorithmic  tangent 
retains  the  characteristic  quadratic  rate  of  convergence  typical  of  a  Newton  method 
when  sufficiently  close  to  the  solution.  Note  that  for  the  last  two  iterations  the 
elements  of  the  discretization  parameters  become  close  to  zero.  Thus,  as  shown 
in  [34],  the  difference  approximation  to  the  algorithmic  tangent  satisfies  the  property 


(5.62)  and  displays  the  essential  characteristics  of  the  Newton  method. 

Next,  the  algorithm  is  tested  for  convergence  under  a  shear  stress  field. 
However,  the  Newton  method  presented  in  Box  3  is  modified  to  streamline 
the  efficiency  of  the  macroscopic  level.  Recognizing  that  the  majority  of  the 
computational  expense  lies  in  each  call  to  the  microscopic  level,  the  previous 
time  step’s  converged  tangent  is  used  for  the  first  two  iterations  of  the  new  time 
step  to  establish  the  auxiliary  points  and  A^.  Also,  the  tangent  operator 
is  updated  every  other  iteration  instead  of  at  every  iteration.  When  using  the 
approximation  (5.65),  this  modification  utilizes  a  minimum  of  three  fewer  calls  to 
the  micromechanical  level  as  outlined  in  Box  3  for  every  two  iterations.  When 
employing  equation  (5.66),  the  modification  uses  one  fewer  call  for  every  iteration. 

For  this  convergence  test,  the  final  configuration  resulting  from  the  one-step 
solution  of  Table  5  (i.e.,  under  a  cumulative  normal  stress  of  =  0^22  =  —10.0  * 
Aq  =  —151.9936)  is  used  as  the  initial  configuration.  From  this  initial  condition,  the 
control  volume  experiences  a  total  shear  stress  of  Ad’j2  =  —3.0  *  Aq  =  —45.5981 
applied  in  two,  three,  four,  six,  12,  and  30  steps.  Table  6  shows  the  results  of  the 
convergence  study  and  suggests  that  the  algorithm  is  convergent  in  the  sense  that 
there  exists  a  macroscopic  shear  motion  to  which  the  solution  tends  as  the  number 
of  steps  increases.  Iterations  for  the  two-step  and  six-step  solutions  failed  to  provide 
a  convergent  solution  because  the  micromechanical  level  did  not  find  an  equilibrium 
configuration  at  the  start  of  a  time  step.  Table  6  also  gives  the  convergence  profile 
for  the  two-step  solution.  This  profile  again  displays  rapid  convergence  in  the 
neighborhood  of  the  solution  but  with  less  computational  effort  than  the  Newton 
method  in  Box  3. 

Recall  that  the  previous  time  step’s  tangent  operator  is  employed  for  the  first 
iteration  of  the  next  time  step.  In  terms  of  the  discretized  Newton  method  described 
in  Box  3,  this  estimate  provides  the  auxiliary  point  A^  which  in  turn  is  used  with  A® 
to  approximate  the  algorithmic  tangent  operator  according  to  either  (5.65)  or  (5.66). 
If  the  previous  step’s  converged  tangent  operator  produces  a  poor  approximation 
of  the  variation  of  the  overall  motion  with  respect  to  the  macroscopic  stress  for 
the  new  time  step,  the  micromechanical  level  may  not  converge  to  an  equilibrium 
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Table  6.  Convergence  test  for  the  196-particle  assembly:  shear 
{k}\{  =  kx  =  I  y.  lO'^units). 


Cumulative  incremental  stress,  Acri2  =  — 1.5  *  Aq  =  —22.7990  : 


Stress  increment 

AAi2,  % 

-1.50 

-0.26453 

-1.00 

n/a 

-0.75 

-0.26434 

-0.50 

-0.26429 

-0.25 

-0.26425 

-0.10 

-0.26423 

Cumulative  incremental  stress,  A<7i2  =  —3.0  *  .Aq  =  —45.5981  : 

Stress  increment 

AAi2,  % 

-1.50 

n/c 

-1.00 

-0.62574 

-0.75 

-0.62614 

-0.50 

n/c 

-0.25 

-0.62698 

-0.10 

-0.62727 

n/a  =  not  applicable;  n/c 

=  no  convergence 

Convergence  profile  for  stress  increment,  A<ti2  =  — 1.5  *  =  —22.7990  : 

Iteration 

Residual 

1 

1.5000000 

2 

9.6489208  x  lO-^ 

3 

1.2397293  x  lO'^ 

4 

4.0669135  X  10"^ 

5 

1.5293010  X  10"“ 

6 

1.1355753  X  10-“ 

configuration.  Figure  27  shows  the  highly  non-linear  relationship  between  the  shear 
stress  ^22  and  the  overall  motion  An  for  the  smallest  stress  increment  of  the  shear 
convergence  test.  This  plot  also  contains  the  data  from  the  non-converged  six-step 
experiment.  The  initial  estimate  of  the  tangent  operator  yields  a  uniform  motion 
which  when  combined  with  the  converged  particle  configuration  of  the  previous 
time  step  causes  the  micromechanical  level  not  to  converge.  Thus,  it  is  important 


to  reiterate  that  whether  the  exact  tangent  operator  or  some  approximation  is 
employed,  the  inverse  algorithm  is  at  best  as  stable  as  the  micromechanical  level 
strain  driven  algorithm. 

Using  the  same  initial  configuration  and  tangent  updating  modification  as  in 
the  shear  convergence  test,  the  particle  assembly  now  experiences  a  complex  stress 
history  consisting  of  alternating  shear  stress  Aa’^2  —  —0.10  *  Aq  =  —1.5199 

and  isotropic  compression.  Figure  28  shows  the  plot  of  the  overall  shear  stress 
d'j2  versus  the  macroscopic  shear  motion  A12.  When  the  micromechanical  level  no 
longer  converges,  the  control  volume  experiences  two  steps  of  isotropic  compression 
Ao’j'j  =  Act22  =  —5.0  *  Ao  =  —75.9968.  Figure  29  shows  the  variation  of 

the  localization  function  minima  with  the  shear  stress.  Despite  the  stress  point 
reaching  a  yield  plateau  during  the  second  shear  loading,  the  localization  function 
does  not  become  zero  or  negative  during  this  experiment.  The  deformed  control 
volume  configuration  at  the  end  of  this  test  appears  in  Figure  30. 
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Figure  21.  Initial  configuration  for  the  64-particle 
closest  packing  control  volume;  particle  radius  = 
1.0  units. 
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NORMALIZED  SHEAR  STREI 


Figure  22.  Normalized  overall  shear  stress  al^lkr 
versus  overall  shear  motion  A\2  showing  formation 
of  a  yield  plateau. 
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Figure  23.  Deformed  configuration  of  the  64- 
particle  closest  packed  control  volume  at  the  end 
of  the  shear  loading  of  Ad'l2  =  —10.0  *  Aq  = 

—151.9936;  particle  radius  =  1.0  units. 
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Figure  24.  Normalized  overall  shear  stress 
versus  overall  shear  motion  A12  showing  the 
effects  of  the  moduli  ratio  0  =  H'  Ikx- 
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1.20E+11 


Figure  26.  Localization  function  versus  orienta¬ 
tion  of  the  plane  of  discontinuity,  at  the  end  of 
the  first  loading  in  the  cyclic  test. 
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Figure  28.  Overall  shear  stress  ^22  versus  overall 
shear  motion  A12  for  shear  loading  followed  by 
isotropic  loading. 
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Figure  29.  Localization  function  minimum  versus 
overall  shear  stress  al2- 
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Figure  30.  Deformed  configuration  for  the 
196-particle  assembly;  mean  particle  radius  = 
0.0092  units. 
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Chapter  6 

Elastic  Contact  Laws  and  Validation 


6.1.  Introduction 

In  1882,  Heinrich  Hertz  introduced  the  basic  premise  used  in  contact  stress 
theory  —  two  bodies  in  contact  can  each  be  regarded  as  an  elastic  half  space  in  the 
proximity  of  the  contact  region.  The  highly  concentrated  contact  stresses  can  then 
be  treated  separately  from  the  general  distribution  of  existing  stresses  which  arise 
from  the  shape  of  the  bodies  and  the  manner  in  which  they  are  supported.  Hertz’s 
contact  theory  relies  upon  a  set  of  basic  assumptions  which  may  be  summarized  as 
follows: 

i.  the  surfaces  are  continuous  and  non-conforming:  the  significant  dimension 
of  the  contact  area  a  is  much  less  than  the  the  relative  radius  of  curvature 
R  where  1/R  =  l/i?i  +  I/R2  and  Ri  and  R2  represent  the  radii  of  the 
bodies  in  contact  (i.e.,  a  ^  R)] 

ii.  the  strains  are  small; 

Hi.  each  solid  is  considered  as  an  elastic  half-space:  a  <C  ,  a  R2  ,  a  <^l 
where  /  denotes  the  dimension  of  the  bodies  in  depth; 

iv.  the  surfaces  are  frictionless. 

Assumption  (f)  ensures  that  the  surfaces  just  beyond  the  contact  region  approximate 
a  plane  surface  of  a  half-space.  Assumption  (u)  guarantees  that  the  strains  within 
the  contact  region  remain  within  the  realm  of  the  theory  of  elasticity.  Assumption 
(fn')  ensures  that  the  stress  field  calculated  on  the  basis  of  a  infinite  solid  will  not 
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be  significantly  influenced  by  the  location  of  its  boundaries  to  the  highly  stressed 
region.  Finally,  assumption  [iv)  is  made  to  guarantee  that  only  a  normal  pressure 
is  transmitted  between  the  contacting  bodies. 

Hertz’s  contact  theory  will  be  employed  to  derive  contact  laws  for  the 
normal  and  tangential  interaction  between  two  infinitely  long  cylinders.  A  secant 
approximation  will  be  used  to  transform  Hertz ’s  nonlinear  contact  expressions  into  a 
linear  contact  law  dependent  upon  the  material  properties  of  the  contacting  bodies. 
The  resulting  contact  stiffnesses  and  kx  will  then  be  substituted  in  the  strain- 
driven  or  stress-driven  algorithm.  The  ability  of  these  algorithms  to  predict  typical 
soil  behavior  will  be  verified  by  subjecting  a  196-particle  control  volume  to  a  specific 
loading  history.  The  results  of  the  numerical  simulation  will  then  be  compared  with 
similar  plane  strain  experimental  data. 


6.2.  Cylindrical  bodies  in  contact:  normal  interaction 

Consider  two  cylindrical  bodies  lying  parallel  to  their  axes  and  pressed  together 
by  a  force  P  as  shown  in  Figure  31.  The  bodies  make  contact  over  a  long  strip  of 
width  2a  parallel  to  the  y-axis.  The  center  of  cylinder  1  and  the  center  of  cylinder  2 
move  towards  the  x  —  y  plane  by  displacements  6i  and  62,  respectively.  If  the 
cylinders  did  not  deform,  their  profiles  would  overlap  as  shown  in  the  insert  of 
Figure  31.  The  displacements  S\  and  62  which  are  measured  positive  into  each 
.cylinder  may  be  found  as  the  limiting  case  of  an  elliptical  contact  area  where  the 
major  axis  of  contact  corresponds  to  the  long  axis  of  the  cylinders  and  becomes  large 
compared  with  the  semi-contact  width  a  as  shown  in  Mindlin  [39].  The  contact  of 
two  long  cylindrical  bodies  then  becomes  two-dimensional  in  nature. 

However,  the  cylinder  contact  problem  may  be  simplified  from  the  outset  by 
employing  assumptions  i  —  iv  in  the  manner  of  Poritsky  [40].  Hertz’s  assumptions 
i  —  iv  neglect  the  curvature  of  the  cylinder’s  boundaries  except  near  the  contact 
region  such  that  the  spreading  of  the  stress  takes  place  in  the  same  manner  as  in  a 
semi-infinite  solid.  The  two-dimensional  plane  strain  contact  configuration  can  then 
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Figure  31.  Cylindrical  bodies  in  contact  under  a  normal  load  P. 


be  transformed  into  an  equivalent  half-space  loaded  by  the  Hertz  contact  pressure 
distribution.  The  distribution  of  pressure  can  be  written  as 

•  p(x)  =  -x‘^)  ,  |x|  <  a;  (6.1) 

where 


SPR 
■kE*  ’ 


R 


Ri^  R2  '' 


1  1  —  vf  1  —  '*^2 

E,  Ei  ' 


(6.2) 

(6.3) 

(6.4) 


#  and  where  P  represents  the  compressive  force  between  the  cylinders  per  unit  axial 

length,  R  is  the  relative  radius  of  curvature,  and  E*  denotes  the  composite  modulus 
of  the  cylinders  in  which  and  En  represent  the  Poisson’s  ratio  and  Young’s 
modulus  respectively  of  cylinder  n. 
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Now,  consider  cylinder  2  in  Figure  31.  For  z  >  0,  the  cylinder  may  be  replaced 
by  a  semi-infinite  solid  loaded  by  the  pressure  distribution  (6.1).  The  semi-infinite 
solid  retains  the  cylinder’s  Poisson’s  ratio  V2  and  Young’s  modulus  E2  ■  Following 
Poritsky  [40],  the  normal  deflection  62{x)  of  the  half  space  along  z  =  0  in  the 
direction  of  the  pressure  can  be  written  as 

62(0:)  = /  \/a2  -s^ln  lx  -s|ds  .  (6.5) 

7rJb2  J-a 


Evaluation  of  this  integral  yields 


Hx)  = 


+  N<a 

+  +<^2,  x>a-, 


where 


w 


=  —  [x  -j-  \/x2  —  a2 
a  L 


(6.6) 


(6.7) 


The  constants  Ci  and  C2  can  be  evaluated  by  considering  the  displacement  62 
relative  to  some  reference  point  where  62  can  be  considered  vanishingly  small. 
Ideally,  the  displacement  would  be  evaluated  relative  to  the  region  at  infinity.  In 
the  case  of  a  finite  object  such  as  a  cylinder,  however,  the  displacement  must  be 
evaluated  relative  to  some  reference  point  in  the  interior  of  the  cylinder.  As  seen 
from  the  logarithmic  term  in  equation  (6.6),  the  displacement  increases  slowly  with 
distance  causing  the  overall  displacement  not  to  be  too  sensitive  of  the  choice  of  the 
reference  point.  Hertz’s  assumptions,  namely  that  a  *C  R2-,  guide  the  choice  of  the 
reference  point  as  equal  to  the  radius  of  cylinder  2.  The  constants  C\  and  C2  then 
become 


Cl  =  -  In 


2R2 


C2  = 


4P(1  -u|) 
TtE2 


1 

a  2’ 

2^  r  2R2 


In 


a 


(6.8) 

(6.9) 


Now,  the  shortening  of  the  branch  vector  connecting  the  centers  of  cylinders  1  and 
2  can  be  evaluated.  As  seen  from  the  particle  kinematics  developed  in  Chapter  2, 
the  compression  of  the  contact  branch  vector  is  of  special  interest  for  developing  the 


particle  contact  constitutive  relationship.  The  shortening  of  the  branch  vector  lying 
within  cylinder  2  can  be  evaluated  from  equation  (6.6)  to  yield: 


^2(0)  = 


2\  r 


TrE2 


1 

a  Z 


(6.10) 


A  similar  expression  can  be  obtained  for  the  compression  of  the  branch  vector  within 
cylinder  1. 

The  total  compression  6  of  the  two  cylinders  equals  the  summation  of  their 
individual  displacements  6i  and  62  as  seen  in  Figure  31.  The  total  compression  of 
the  branch  vector  can  be  written  as 


<5(0)  =  (5i(0)  + (52(0) 
_  4P(1  -u^)  r 
•kE 


(6.11) 

(6.12) 


a  a 

This  expression  for  the  relative  approach  of  the  centers  of  two  infinitely  long  cylinders 
represents  a  nonlinear  relation  in  terms  of  the  contact  force  per  unit  length  P.  The 
nonlinearity  arises  from  the  dependence  of  the  semi-contact  width  a  on  the  contact 
force  P  as  seen  in  equation  (6.2). 


6.3.  Cylindrical  bodies  in  contact:  tangential  interaction 

When  considering  a  tangentially  loaded  point  of  contact,  a  distinction  must 
be  made  between  slip  and  sliding.  Slip  occurs  when  there  is  a  relative  tangential 
displacement  over  only  a  portion  of  the  contact  surface.  The  region  within  which  no 
relative  tangential  displacement  occurs  is  termed  the  istick’  region.  In  contrast, 
sliding  occurs  when  tangential  force  reaches  a  critical  value  to  cause  relative 
displacement  over  the  whole  contact  surface.  The  particular  case  considered  here 
will  be  elastic  bodies  in  contact  under  the  condition  of  no  slip.  When  the  tangential 
traction  does  not  exceed  the  limiting  value,  the  no  slip  condition  requires  that  all 
surface  points  within  the  ‘stick’  region  undergo  the  same  tangential  displacement. 

The  specific  case  of  interest  consists  of  the  compression  of  two  infinitely  long 
cylinders  by  a  normal  force  P  per  axial  length  and  the  subsequent  application  of  a 
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Figure  32.  Tangential  displacements  of  two  contacting  cylinders. 

tangential  force  T  per  unit  length.  The  Hertz  contact  theory  of  the  previous  section 
gives  the  pressure  distribution  arising  from  the  force  P,  the  semi-contact  width  a, 
and  the  total  normal  displacement  6  of  the  cylinders.  The  complete  contact  width 
—a<x<a  represents  the  ‘stick’  area  where  no  slip  occurs.  Figure  32  shows 
the  tangential  force  T  at  the  contact  interface  and  the  resulting  tangential,  elastic 
displacements  of  the  contact  strip.  The  tangential  displacements  71  and  71  represent 
the  tangential  displacement  of  the  center  of  cylinder  1  and  2  respectively.  Note  that 
the  center  of  each  cylinder  has  been  considered  distant  from  the  loaded  region. 

Mindlin  [39]  recognized  that  if  no  slip  occurs  at  the  contact  surface,  no  change 
in  the  normal  component  of  traction  occurs  across  the  surface.  Thus,  when  the 
bodies  have  the  same  elastic  properties,  the  shape  and  size  of  the  contact  area 
and  the  distribution  of  the  normal  pressure  will  be  independent  of  the  tangential 
force  [41].  The  profiles  of  the  two  contacting  surfaces  fix  the  shape  and  size  of  the 
contact  area.  The  stresses  and  deformations  resulting  from  a  normal  pressure  and 
the  tangential  traction  will  be  independent  of  each  other,  and  superposition  can  be 
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employed  to  find  the  resultant  stress. 

Consider  two  cylinders  1  and  2  in  contact  over  a  region  |a:|  <  a  on  2:  =  0  and 
subject  to  a  tangential  traction  given  by 

t{x)  =  ~  ’  l®l  <  (6.13) 

Again,  the  assumptions  of  Hertz  can  be  employed  to  transform  the  problem  into  an 
elastic  half  space  loaded  by  the  tangential  traction  (6.13)  over  |x|  <  a. 

The  tangential  displacements  71  and  72  can  be  evaluated  by  employing  their 
basic  analogy  with  the  normal  displacements  61  and  82-  Under  the  action  of  identical 
distributions  of  tangential  and  normal  tractions  t{x)  and  p(a:),  the  tangential 
displacements  71  and  72  analogous  with  the  normal  displacements  and  82-  For 
example,  the  normal  tractions  given  in  (6.1)  along  the  plane  boundary  of  the 
half  space  can  be  replaced  by  the  tangential  tractions  (6.13),  and  the  normal 
displacements  and  82  found  in  the  previous  section  can  be  replaced  by  the 
tangential  displacements  71  and  72  respectively.  The  analogy  assumes  that  the 
same  reference  point  has  been  used  to  evaluate  the  constants. 

By  employing  the  normal-tangential  analogy,  the  solution  for  the  relative 
approach  of  the  centers  of  the  cylinders  given  in  equation  (6.12)  can  be  employed 
to  write  the  total  tangential  displacement  7  of  the  centers  of  the  cylinders  as 


7(0)  =  71(0)  -72(0) 
4T(1  -t;2) 


ttE 


(6.14) 

(6.15) 


Note  that  since  the  cylinders  experience  mutually  equal  and  opposite  tangential 
tractions,  their  displacements  72(6)  and  71  (0)  in  (6.14)  will  be  of  opposite  signs. 

Mindlin  [39]  considered  the  general  case  of  a  normal  and  tangential  forces 
applied  across  an  elliptic  contact  surface  of  a  pair  of  elastic  bodies.  For  example, 
the  solution  of  the  tangential  displacements  for  an  elliptical  contact  assuming  no 
slip  is  given  by 


7  = 


ttG 


K-^(K-E) 


(6.16) 


and  G  denotes  the  shear  modulus,  K  and  E  are  complete  elliptic  integrals  of  the  # 

first  and  second  kind,  respectively,  and  £  is  the  axial  length  of  the  contact  area.  The 

contact  stresses  for  the  two-dimensional  case  of  contact  between  two  cylinders  may 

be  obtained  from  these  solutions  as  a  limiting  case  by  allowing  the  contact  area  to 

become  infinitely  long  along  the  long  axis  of  the  cylinders.  The  resulting  relations  • 

for  the  normal  and  tangential  displacements  are  then  given  by  equations  (6.12)  and 

(6.15). 


6.4.  Linearization  of  the  contact  laws 

The  contact  laws  for  the  normal  and  tangential  interaction  of  long  cylindrical 
bodies  given  by  equations  (6.12)  and  (6.15)  represent  nonlinear  relationships  between 
the  contact  force  and  displacement.  The  nonlinearity  arises  from  the  displacement’s 
logarithmic  variation  with  the  normal  contact  force  P.  Simplified  expressions  for 
the  normal  and  tangential  contact  laws  can  be  derived  with  a  secant  approximation. 

Consider  minimum  and  maximum  normal  contact  forces  Pmin  t^nd  Pmax-  The 
slope  of  the  line  connecting  these  two  points  on  the  nonlinear  contact  laws  given  in 
equations  (6.12)  and  (6.15)  defines  the  secant  approximation.  For  choosing  a  value 
of  the  minimum  normal  contact  force,  the  condition  of  incipient  contact  requires 
Pfnin  =  0.  Two  criterion  guide  the  choice  of  the  maximum  normal  contact  force 
Ptnax-  First,  the  Hertz  assumptions  i—iv  demand  that  Pmax  satisfy  the  requi^ements 
of  small  strain,  namely  a  ^  R.  Second,  the  accuracy  of  the  secant  approximation 
demands  that  Pmax  he  similar  to  the  actual  maximum  normal  force  between  the 
two  cylindrical  bodies.  With  Pmin  =  0  and  the  appropriate  choice  for  Pmax  that 
satisfies  these  criterion,  the  secant  approximation  to  the  Hertzian  contact  laws  given 
in  equations  (6.12)  and  (6.15)  take  the  form 

P  =  (6.18) 

T  =  kri  (6.19) 

where  and  fcj'  denote  the  slopes  of  the  normal  and  tangential  linear  relations, 
respectively.  Physically,  kj^  and  Arj  may  be  interpreted  as  linear  normal  and 
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tangential  contact  spring  constants. 

Consider  expressions  for  the  normal  and  tangential  contact  stiffnesses  kff  and 
kT-  The  normal  Hertzian  contact  law  (6.12)  repeated  here  for  convenience  as 


4(1  —  t!^) 

ttE 


In 


i?l  R2T^E* 
2RP 


P 


(6.20) 


where  the  semi-contact  width  a  (6.2)  has  been  employed.  Substituting  the  minimum 
and  maximum  normal  forces  {0,Pmax)  iiito  this  expression  defines  the  slope  of  the 
normal  contact  approximation  and  gives 


kN  = 


'kE 

4(1  -u2) 


In 


R\R2T^  E* 

^RPmax 


(6.21) 


Thus,  the  normal  secant  approximation  (6.18)  can  be  written  as 


P  = 


irE 


4(l-t;2)[  2RP, 


RiR2'rcE* 

In  .  „  - - h  1 
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-1 


fcjv 


(6.22) 


In  a  similar  manner,  the  secant  approximation  to  the  Hertzian  tangential  law  can 
be  written  as 


r  = 


irE 

4(1  -i;2) 


In 


R\R2^  E* 
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-1 
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(6.23) 
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Note  that  the  normal  and  tangential  contact  stiffnesses  are  equal  and  depend  upon 
the  material  properties  and  radii  of  the  contacting  cylinders. 


Figure  33  contains  a  comparison  of  the  nonlinear  Hertz  contact  law  and  its 
secant  approximation  (6.22)  for  similar  cylinders.  The  cylinders  possess  identical 
material  and  geometric  properties.  The  basis  of  the  comparison  consists  of  the 
nondimensionalized  contact  force  Pl{kR\)  and  the  normal  indentation  SjRi  where 
k  =  7rFl/4(l  —  v^)  and  R\  =  R2  =  1.0  consistent  units  implied.  The  range  of  the 
secant  approximation  has  been  restricted  to  the  region  where  0  <  Pjk  <  1/625. 
In  general,  over  the  region  of  interest  the  linearized  contact  law  provides  a  good 
approximation  to  the  Hertzian  law. 


Expressions  for  the  contact  forces  Tn-v\  discussed  in  Chapter  2  can  be  obtained 
from  the  discretization  of  the  linearized  contact  laws  (6.22)-(6.23).  Assume  that  at 
time  in  the  contacting  cylinders  A  and  B  —  particles  in  plane  strain  —  have  well 
defined  centroids  a;^  and  .  The  contact  force  P  between  two  cylinders  can  now 
represent  the  normal  component  of  the  contact  force  ///at  time  tn+i-  The  normal 
contact  law  can  now  be  written  in  incremental  form  as 

fn+l  =  -He{Sn+l)kN6n+l  (6.24) 

where  the  ramp  function  Hf(Sn+i)  is  as  defined  in  Chapter  2,  the  overlap 
equals  the  shortening  of  the  branch  vector  <^(0),  and  fc/vr  denotes  the  normal  contact 
stiffness  as  defined  by  (6.21). 

For  the  tangential  interaction,  the  contact  stiffness  in  (6.23)  defines  the  elastic 
tangential  force.  When  the  contact  mode  is  ‘stick,’  the  tangential  component  of  the 
contact  force  fx  at  time  tn+i  can  be  written  as 

(/r)n+l  =  (/T)n  +  H,{S„+i)kTA7  (6.25) 

where  the  A7  represents  the  increment  of  tangential  displacement.  Note  that 
the  tangential  displacement  7  of  (6.23)  represents  the  total  elastic  tangential 
displacement  at  the  contact.  If  the  contact  is  in  ‘sliding’  mode,  the  tangential 
force  then  becomes 


(/r)«+i  =  {fT)n  +  ife(<^«+i)fcr(A7  -  A7^)  (6.26) 

The  remaining  results  of  Chapter  2  apply  for  these  contact  modes  assuming  that 
over  the  finite  incremental  motion  the  particles  remain  in  contact  throughout  the 
time  step.  Now,  the  incremental  normal  and  tangential  contact  laws  (6.24)-(6.26) 
form  the  basis  of  the  contact  model  between  two  particles.  Thus,  the  Hertz  solution 
to  the  contact  of  infinite  circular  cylinders  forms  the  foundation  of  the  contact  model 
underlying  the  numerical  algorithms  of  Chapter  4  and  Chapter  5. 


6.5.  Numerical  simulations 


With  the  addition  of  the  Hertz  solution  to  the  contact  model,  the  behavior 
of  the  numerical  material  takes  on  the  added  physical  meaning  of  the  response  of 
a  similarly  graded  granular  material  subject  to  plane  strain  loading.  Herein,  the 
circular  cylinders  represent  the  action  of  individual  soil  particles  contained  in  a 
repeating  control  volume  subject  to  plane  strain  loading.  The  predicted  behavior 
of  the  numerical  material  can  now  be  compared  with  experimental  plane  strain  test 
data. 

Before  comparing  the  numerical  model  and  experimental  tests,  one  must 
consider  and  clarify  the  limitations  and  shortcomings  of  such  a  comparison.  Consider 
the  plane  strain  basis  of  the  comparison.  The  numerical  model  represents  a  classical 
plane  strain  formulation  with  infinitely  long  cylinders  identified  as  soil  particles. 
However,  actual  three-dimensional  tests  model  plane  strain  behavior  by  restricting 
deformation  in  the  longest  specimen  direction.  The  typically  tabular  shaped  samples 
contain  a  very  large  number  of  finite  sized  particles.  Thus,  the  comparison  reduces 
to  a  purely  plane  strain  numerical  model  and  a  three-dimensional  experiment 
configured  to  model  plane  strain  behavior. 

Also,  consider  the  boundary  condition  basis  of  the  comparison.  Recall  that 
the  assumption  of  homogeneous  deformation  is  implicit  in  the  numerical  model. 
The  basis  for  comparison  of  the  numerical  model  and  experimental  results  rests  on 
the  experimental  apparatus’  ability  to  allow  homogeneous  deformation.  However, 
such  effects  such  as  end  plate  friction  and  membrane-soil  interaction  are  inherent 
in  physical  experiments.  The  prevalence  of  these  factors  in  experimental  results 
directly  influences  the  applicability  of  the  homogeneous  deformation  assumption  for 
the  experiment. 

Finally,  consider  the  possible  scaling  effects  produced  by  using  a  relatively  small 
number  of  particles  to  represent  a  soil  sample.  The  numerical  model  assumes  that 
the  particles  in  the  control  volume  are  sufficiently  representative  to  capture  the 
geometric  characteristics  of  the  soil.  Yet,  a  large  number  of  particles  are  typically 
necessary  to  capture  a  soil’s  gradation,  packing,  and  overall  void  structure.  Thus, 


a  trade  off  typically  exists  between  the  cost  of  the  numerical  simulation,  i.e.,  the 
number  of  particles  in  a  control  volume,  and  how  representative  the  control  volume 
will  be  of  the  model  soil. 

Experimental  results 

Many  investigators  [42-47]  have  studied  the  plane  strain  behavior  of  granular 
materials.  Vardoulakis  [48]  and  Vardoulakis  and  Graf  [49]  present  extensive  studies 
of  the  plane  strain  compression  of  dry  sand  that  appear  representative  of  the 
behavior  of  granular  materials  in  plane  strain  loading.  Their  test  samples  appear 
tabular  in  shape  with  their  long  axis  perpendicular  to  the  loading  direction.  The 
sample  is  placed  in  a  conventional  triaxial  loading  frame  with  two  stiff  frictionless 
plates  to  restrict  the  strain  to  two  dimensions.  The  cell  pressure  remains  constant 
while  the  loading  piston  compresses  the  sample  under  constant  driving  speed.  The 
piston  displacement  and  force  are  measured  throughout  the  test.  Vardoulakis 
and  Goldscheider  [50]  explain  in  greater  detail  the  biaxial  apparatus  and  testing 
procedure. 

Two  different  experimental  results  will  be  used  for  comparison  with  the 
numerical  model.  The  first  test  from  Vardoulakis  [48]  consists  of  the  plane  strain 
compression  of  a  dry  Osterchelde  sand  at  a  confining  pressure  of  29.43  N/cm  .  The 
sample’s  initial  void  ratio  equaled  e  =  0.621  and  specific  gravity  was  Gs  =  2.66.  The 
second  test  from  Vardoulakis  and  Graf  [49]  was  performed  at  a  confining  pressure 
of  19.62  N/cm^  on  a  dry  sample  of  Karlsruhe  sand  with  an  initial  void  ratio  equal 
,to  e  =  0.577. 

The  reported  experimental  data  for  both  tests  consists  of  the  variation  of  the 
axial  force  Pa  with  the  axial  displacement  Ua-  For  comparison  with  the  numerical 
simulations,  the  deviator  stress  Actq  and  the  axial  strain  Ca  have  been  calculated 
with  the  following  equations: 

A<r.  =  a,  -  as  =  (6.27) 

^  (6-28) 
where  Aq  represents  the  initial  unloaded  cross-sectional  area  of  the  sample,  and  la 


denotes  the  sample’s  initial  length  in  the  direction  of  loading.  The  variation  of  the 
deviator  stress  Acto  and  axial  strain  Cq  can  now  be  plotted. 

Numerical  material  properties 

The  two  control  volumes  used  in  the  numerical  simulations  consist  of 
parallelepipeds  of  infinite  circular  cylinders.  Each  volume  contains  49  and  196 
particles  (cylinders)  randomly  arranged.  All  loading  takes  place  perpendicular  to 
the  axis  of  the  cylinders.  Recall  that  the  control  volumes  represent  one  cell  within 
a  periodic  structure.  The  control  volume  of  interest  is  surrounded  in  all  directions 
by  identical  parallelepipeds.  All  particles  within  the  cell  possess  the  same  material 
properties  but  not  necessarily  the  same  geometric  properties.  The  particles  may  be 
randomly  dispersed  throughout  the  cell  and  have  varying  radii. 

In  choosing  the  material  properties,  namely,  the  Young’s  modulus  and  Poisson’s 
ratio,  of  the  repeating  cell  particles  consideration  must  be  given  to  the  composition 
and  origin  of  soil  particles.  The  principal  constituents  of  the  solid  phase  of  soils  are 
various  amounts  of  crystalline  and  noncrystalline  clay  materials,  nonclay  minerals, 
and  precipitated  salts.  In  particular,  nonclay  materials  such  as  sands,  gravels,  and 
the  majority  of  silts  are  typically  comprised  of  rock  fragments  or  mineral  grains  of 
the  common  rock-forming  minerals.  Thus,  the  properties  of  the  pre-existing  rock 
from  which  the  individual  grains  originate  can  be  used  to  estimate  the  material 
properties  of  the  particles  in  the  numerical  model. 

Consider  numerical  soil  particles  that  originate  from  rock  fragments  or  mineral 
grains  of  common  sedimentary  rocks  such  as  sandstone,  siltstone,  shale,  and 
claystone.  The  particle  size  of  the  reference  materials  varies  from  .2  —  .00625  cm  for 
sandstone,  .00625  — .00039  cm  for  siltstone,  and  <  .00039  cm  for  shale  and  claystone. 
The  Young’s  modulus  and  Poisson’s  ratio  of  these  common  quartz  materials  have 
been  determined  from  unconfined  compression  tests.  The  rock  samples  are  typically 
core  specimens  of  length  approximately  twice  the  diameter.  A  summary  of  these 
properties  for  the  reference  materials  and  many  other  common  rock  materials  can 
be  found  in  Physical  Properties  of  Rocks  and  Minerals  [51].  The  Young’s  modulus 
and  Poisson’s  ratio  of  the  numerical  soil  particles  are  chosen  to  lie  within  the  range 


Table  7.  Model  and  material  parameters. 


Model  parameters 

Macromechanical  error  tolerance;  RTOL  =  1.0  x  10“^. 
Micromechanical  error  tolerance;  rtol  =  1.0  x  10“®. 
Material  constants 

Particle  Young’s  modulus;  E  =  0.25  —  2.0  GPa. 
Particle  Poisson’s  ratio:  v  =  0.20. 

Particle  friction  angle:  (j>  =  30°. 

Particle  contact  cohesion;  ao  =  0. 

Contact  hardening  parameter:  H'  =  0. 


of  the  reference  material’s  properties. 

Table  7  summarizes  the  model  and  material  parameters  of  interest.  The 
error  tolerances  listed  represent  the  convergence  criterion  used  in  the  numerical 
algorithms  described  in  Chapter  4  and  Chapter  5.  A  Young’s  modulus  in  the 
range  of  0.25  —  2.0  GPa  and  a  Poison’s  ratio  equal  to  0.20  has  been  used.  All 
numerical  simulations  assume  an  elastic-perfectly  plastic  contact  model  for  each 
particle  contact  by  employing  a  contact  hardening  parameter  equal  to  if'  =  0. 

Forty  nine-particle  control  volume 

The  initial  control  volume  configuration  for  this  example  is  shown  in  Figure  34. 
The  control  volume  configuration  isolates  the  effect  of  a  random  arrangement  of 
particles  from  the  effect  of  randomly  sized  particles.  All  particles  in  the  repeating 
cell  have  a  radius  equal  to  0.0092  cm.  This  particular  radius  has  been  chosen  to 
equal  the  mean  particle  radius  of  a  randomly  sized  and  randomly  arranged  196- 
particle  control  volume.  The  voids  in  the  particle  configuration  seen  in  Figure  34 


have  been  created  by  removing  seven  particles  from  a  closest  packed  arrangement 
of  56  particles.  The  initial  cell  contains  126  particle  contacts  and  has  dimensions  of 
0.12880  cm  x  0.12748  cm. 


For  all  simulations,  the  value  of  the  maximum  normal  contact  force  between 
two  particles  equals  Pmax  =  2.0  N  for  all  particle  contacts  in  the  repeating  cell. 
Thus,  the  normal  and  tangential  contact  stiffnesses  remain  constant  for  all  particle 
contacts.  Also,  the  ramp  function  e  equals  zero.  Recall  that  e  defines  a  finite 
transition  zone  from  no-contact  to  full-contact.  With  e  =  0,  an  exact  simulation  of 
the  abrupt  transition  to  full  contact  can  be  obtained.  As  a  result,  the  ramp  function 
He{8n+i)  exactly  models  the  Heaviside  function  and  has  the  form 

Finally,  the  overall  stress  state  of  the  control  volume  at  the  end  of  each  timestep 
has  been  calculated  using  a  volume  equal  to  the  current  cross  sectional  area  times  a 
unit  thickness.  The  stress  field  then  becomes  a  Cauchy  representation  of  the  stress 
state  since  the  current  deformed  configuration  has  been  used. 

2 

The  first  numerical  experiment  reproduces  the  confining  stress  of  29.43  N/cm 
in  the  Vardoulakis  [48]  test  on  dry  Osterchelde  sand.  The  control  volume  experiences 
two  increments  of  isotropic  stress  of  =  Au'22  =  —14.715  N/cm  ;  A5’j2  =  0.0. 
Following  the  isotropic  compression,  the  control  volume  experiences  increments  of 
axial  loading  of  Ad-Jj  =  —0.50  N/cm  ;  A0-22  =  Ao’j2  =  0.0  until  the  algorithm  no 
longer  converges. 

Figure  35  presents  the  response  of  the  control  volume  by  plotting  the  variation 
of  the  deviator  stress  Aa*  =  Ad'j^  —  Ad-22  versus  the  axial  strain  ei  along  with 
the  experimental  results  obtained  by  Vardoulakis  [48].  The  numerical  simulations 
have  been  run  for  four  different  Young’s  moduli  of  0.25,  0.50,  1.00,  and  2.00  GPa. 
Figure  35  indicates  that  as  the  Young’s  modulus  increases  the  response  of  the 
control  volume  becomes  stiffer  as  indicated  by  a  steeper  deviator  stress-axial  strain 
curve.  Figure  36  shows  an  enlargement  of  the  initial  stress-strain  behavior  of  the 
numerical  and  experimental  results.  All  four  values  of  the  Young’s  modulus  provide 
good  correlations  with  the  experimental  results.  Figure  37  shows  the  deformed 
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configuration  of  the  control  volume  at  the  end  of  the  numerical  simulation  with  the 
Young’s  modulus  equal  to  =  0.25  GPa.  At  the  end  of  this  loading,  the  control 
volume  contains  119  particle  contacts  of  which  35  have  plastified,  i.e.,  entered  “slip” 
mode,  and  a  maximum  normal  contact  force  equal  to  3.09856  N. 

Figure  38  presents  the  behavior  of  the  control  volume  to  a  loading  history 
similar  to  the  one  employed  by  the  Vardoulakis  and  Graf  [49]  on  Karlsruhe  sand. 
The  control  volume  experiences  an  isotropic  stress  of  =  Ad'22  —  —9.81  N/ cm  ; 
Ad-J2  =  0.0  applied  in  two  increments  to  provide  a  confining  stress  of  19.62  N/ cm  . 
Next,  the  loading  consists  of  axial  stress  increments  of  Ad-ji  =  —0.50  N/cm^; 
Aa22  =  Aai2  =  0.0.  The  predicted  behavior  at  the  same  four  different  Young’s 
moduli  appear  along  side  the  experimental  test  data.  Like  the  previous  simulation, 
the  predicted  stress-strain  behavior  appears  similar  to  the  experimental  data.  In 
particular.  Figure  39  shows  the  comparison  of  predicted  and  experimental  data  at 
small  strains.  The  deformed  control  volume  at  the  end  of  the  loading  with  Young’s 
modulus  equal  to  E  =  0.25  GPa  can  be  seen  in  Figure  40.  At  this  point,  the 
repeating  cell  contains  119  contact  of  which  37  have  plastified  and  a  maximum 
normal  contact  force  equal  to  2.11516  N. 

One  hundred  ninety  six-particle  control  volume 

The  initial  control  volume  configuration  for  this  example  is  shown  in  Figure  41. 
The  repeating  cell  contains  196  randomly  sized  randomly  arranged  particles 
originally  configured  by  Kuhn  [36].  The  mean  particle  radius  equals  0.0092  cm,  and 
the  initial  cell  dimensions  are  0.2566  cm  x  0.2565  cm.  Unless  otherwise  noted  the 
model  and  material  parameters  equal  those  listed  in  Table  7.  Three  analyses  have 
been  performed  to  evaluate  the  performance  of  the  numerical  model  on  a  randomly 
configured  control  volume.  The  first  analysis  considers  the  sensitivity  of  the  model 
to  certain  parameters  such  as  the  ramp  function  e  and  maximum  normal  contact 
force  Pmax-  The  second  analysis  evaluates  the  convergence  characteristics  of  the 
model.  The  final  analysis  compares  the  deviator  stress-axial  strain  behavior  with 
experimental  results. 


Sensitivity  analysis 

The  numerical  model  performance  is  tested  for  sensitivity  to  variations  in  the 
ramp  function  e  and  maximum  normal  contact  force  Pmax-  For  all  sensitivity 
simulations,  the  control  volume  experiences  a  loading  history  consisting  of  two 
increments  of  isotropic  compression  followed  by  axial  loading  increments  of  = 
—0.50  N/cm^;  Ad’22  =  ^^12  —  0-0- 

In  the  ramp  function  study,  the  two  isotropic  compression  increments  equal 

o 

AdJj  =  Ad'22  =  —14.715  N/cm  ;  Ad-j2  =  0.0,  and  the  Young’s  modulus  and  Pmax 
equal  2.0  GPa  and  2.0  N,  respectively.  Three  different  values  of  e  have  been  defined: 
0,  0.05  *Rmini  and  0.1  *Rmin  where  Rmin  equals  the  radius  of  the  smallest  particle 
in  the  control  volume.  For  the  196-particle  control  volume,  Rmin  equals  0.00503  cm. 
Thus,  the  ramp  function  creates  a  constant  layer  of  softer  material  of  thickness  e/2 
around  each  particle  in  the  control  volume.  For  example,  when  e  =  0.05  *  Rmin, 
every  particle  in  the  repeating  cell  has  a  softer  layer  of  thickness  equal  to  5%  of  the 
smallest  particle  radius  or  0.0025  mm. 

With  the  ramp  function  e  equal  to  zero,  the  model  failed  to  converge  for  the 
initial  increments  of  isotropic  compression.  The  random  sizes  of  the  particles  within 
the  repeating  cell  naturally  cause  relatively  large  particles  to  come  into  contact  with 
smaller  particles.  Since  the  contact  stiffnesses  and  kx  depend  on  the  radii  of 
the  contacting  particles,  the  contact  stiffnesses  throughout  the  repeating  cell  vary 
depending  on  which  particles  lie  in  contact.  Thus,  the  different  sized  particles  in  the 
196-particle  control  volume  produces  a  numerically  “stiff”  system.  This  behavior 
can  be  contrasted  with  the  performance  of  the  49-uniformly  sized  particle  control 
volume.  In  this  case,  the  contact  stiffnesses  fcjy  and  kx  are  constant  for  all  contacts  in 
the  control  volume.  A  ramp  function  equal  to  zero  does  not  hinder  the  convergence 
as  shown  in  Figure  35  and  Figure  38. 

The  performance  of  the  model  with  non-zero  ramp  function  can  be  seen  in 
Figure  42.  The  deviator  stress-axial  strain  behavior  for  e  equal  to  0.05  *  Rmin 
and  0.1  *  Rmin  have  been  plotted  against  the  experimental  data  for  Osterchelde 
sand  [48].  The  figure  shows  the  better  convergence  behavior  of  the  model  as  it 
becomes  less  stiff  with  a  larger  ramp  function.  However,  the  stress-strain  behavior 


does  not  differ  significantly  in  the  initial  portion  of  the  loading  history.  Thus,  a 
larger  ramp  function  enhances  the  convergence  of  the  model  without  significantly 
sacrificing  accuracy.  Yet,  care  should  be  exercised  while  choosing  a  value  of  £.  A 
large  of  a  ramp  function  could  contaminate  the  accuracy  of  the  solution  by  creating 
too  large  of  a  soft  band  around  the  particles. 

The  sensitivity  of  the  model  to  the  choice  of  the  maximum  normal  contact 
force  Pmax  is  investigated  next.  The  two  increments  of  isotropic  compression  equal 
AaJi  =  Acr22  =  -9.81  N/cm^;  Actj2  =  O-O-  Figure  43  contains  the  stress-strain 
behavior  with  a  Young’s  modulus  of  0.25  GPa  and  Pmax  equal  to  35.0,  10.0,  and 
2.0  N.  Also,  indicated  in  the  figure  are  the  values  of  the  actual  maximum  normal 
contact  force  P^*^  in  the  control  volume  at  the  end  of  the  simulation.  Figure  43 
shows  a  similar  stress-strain  behavior  for  different  values  of  Pmax  iii  fhe  early 
portions  of  the  simulations.  As  the  value  of  Pmax  approaches  the  value  of  PmL^ 
the  secant  approximation  to  the  Hertzian  contact  laws  becomes  more  accurate. 
Correspondingly,  the  convergence  properties  of  the  model  improve. 

Convergence  analysis 

The  numerical  model  is  next  tested  for  convergence.  The  Young’s  modulus 
and  ramp  function  for  the  two  convergence  analyses  equal  0.50  GPa  and  0.1  *  Rmin-, 
respectively.  In  the  first  convergence  test,  the  control  volume  initially  experiences 
a  confining  stress  of  29.43  N/cm^  applied  in  two  increments  of  Aa^j  =  Act22  = 
—14.715  N/cm^;  A^*2  =  0.0.  The  second  test  uses  a  confining  stress  of  19.62  N/cm^ 
applied  in  two  increments  of  Aa^j  =  Ad'22  =  —9.81  N/ cm^;  Actj2  —  ^Fe  control 
volume  is  then  compressed  axially  to  a  deviator  stress  of  Aa*  =  15  N/cm  in  3,  30, 
150  increments  of  Actjj. 

Table  8  presents  the  results  of  the  convergence  study.  The  results  indicate  that 
the  model  is  convergent  in  the  sense  that  an  axial  strain  exists  to  which  the  solution 
tends  as  the  number  of  steps  increases.  For  the  simulation  with  a  confining  stress  of 
29.43  N/cm^,  the  error  of  the  three-step  solution  is  in  the  order  of  0.5%  relative  to 
the  150-step  solution.  For  the  test  with  a  confining  stress  of  19.62  N/cm  ,  the  model 
failed  to  converge  for  the  largest  increment  of  axial  strain.  As  the  increment  size 


Table  8.  Convergence  test  for  the  196-paxticle  assembly:  £'  =  0.50  GPa. 


2  2 
Confining  stress:  ajj  =  ^22  =  29.43  N/cm  ;  deviator  stress:  Aa*  =  15.0  N/cm 


Axial  stress  increment  (N/cm^) 

AAii,  % 

-5.0 

-0.38375 

-0.5 

-0.38539 

-0.1 

-0.38555 

2 

Confining  stress:  ajj  =  ^22  =  19.62  N/cm 

2 

;  deviator  stress:  Aa*  =  15.0  N/cm  . 

Axial  stress  increment  (N/ciri?) 

AAii,  % 

-5.0 

n/c 

-0.5 

-0.48953 

-0.1 

-0.48982 

n/c  =  no  convergence 


decreased,  the  results  appeared  convergent  in  the  same  manner  as  the  convergence 
test  at  the  higher  confining  stress. 

Simulations  of  experimental  tests 

The  behavior  of  the  196-particle  control  volume  is  compared  to  the 
experimental  data  for  Osterchelde  and  Karlsruhe  sands.  The  loading  history  for 
these  numerical  simulations  is  identical  to  the  corresponding  tests  on  the  49-particle 
repeating  cell.  Again,  the  effect  of  different  Young’s  moduli  has  been  investigated 
with  E  =  0.25,  0.50,  1.00,  and  2.00  GPa.  The  ramp  function,  however,  now  equals 
0.1  *  Rjnin' 

Figure  44  shows  a  comparison  of  the  deviator  stress- axial  strain  behavior  of 
the  repeating  cell  at  a  confining  stress  of  29.43  N/cm  with  that  of  the  Osterchelde 
sand  [48].  The  predicted  behavior  at  all  four  values  of  the  Young’s  modulus  provide 
good  correlations  with  the  experimental  results  particularly  in  the  initial  portion  of 
the  loading  history.  Figure  45  shows  an  enlargement  of  this  region  of  the  stress- 
strain  graph.  The  simulation  with  the  Young’s  modulus  of  2.0  GPa  provides  the 
best  correlation  with  the  experimental  data.  Figure  46  presents  the  deformed  control 
volume  at  the  end  of  the  test  with  E  =  0.25  GPa.  For  this  configuration,  the  control 


Table  9.  Initial  secant  modulus. 


2 

Confining  stress:  =  ^22  ~  29.43  N/cm  . 


Young’s 
Modulus  (GPa) 

Secant  Modulus  (N/cm^) 
49  Particle  Repeating  Cell 

Secant  Modulus  (N/cm^) 

196  Particle  Repeating  cell 

0.25 

33.7338 

29.8073 

0.50 

59.1248 

38.5258 

1.00 

99.3776 

49.7157 

2.00 

158.3656 

63.1377* 

2 

Experimental  data:  initial  secant  modulus  =  76.8263  N/cm  . 

Confining  stress:  —  ^22  =  19.62  N/cm^. 

Young’s 

Secant  Modulus  (N/cm^) 

Secant  Modulus  (N/cm^) 

Modulus  (GPa) 

49  Particle  Repeating  Cell 

196  Particle  Repeating  cell 

0.25 

33.3035 

23.7078 

0.50 

56.6170 

30.6107 

1.00 

90.7327 

39.0597* 

2.00 

150.4227* 

48.7559* 

2 

Experimental  data:  initial  secant  modulus  =  65.6802  N/cm  . 
*  Secant  modulus  based  on  the  final  axial  stress  and  strain. 


volume  contains  428  particle  contacts  of  which  56  have  plastified,  and  the  maximum 
normal  contact  force  in  the  control  volume  equals  2.07579  N. 

Figure  47  shows  the  stress-strain  behavior  of  the  numerical  and  experimental 
tests  at  a  confining  stress  of  19.62  N/cm^.  The  experimental  data  is  for  a  Karlsruhe 
sand  [49].  Again,  the  predicted  behavior  compares  well  with  the  experimental  data. 
The  observed  trend  of  the  simulations  with  lower  Young’s  modulus  providing  better 
convergence  is  apparent  in  this  simulation.  A  smaller  Young’s  modulus  makes  the 
model  less  numerically  “stiff”  in  a  similar  manner  to  the  effect  of  a  larger  ramp 
function.  Figure  48  shows  the  initial  stress-strain  behavior.  Figure  49  shows  the 
deformed  repeating  cell  at  the  end  of  the  test  with  E  =  0.25  GPa.  Here  the  control 
volume  contains  408  contacts  of  which  60  have  plastified,  and  the  maximum  normal 
contact  force  in  the  control  volume  equals  1.69781  N. 

Table  9  contains  a  comparison  between  the  initial  secant  modulus  of  the 


numerical  simulations  and  experimental  data.  The  secant  modulus  has  been 
calculated  for  an  axial  strain  equal  to  0.5%  for  the  49  and  196-particle  control 
volumes  at  confining  stresses  of  29.43  N/cm  and  19.62  N/cm  .  The  results  indicate 
that  for  the  range  of  the  Young’s  moduli  considered  the  model  produces  an  initial 
secant  modulus  of  the  same  order  of  magnitude  as  the  experimental  data.  In  general, 
the  196-particle  control  volume  performs  better  than  the  more  uniform  49-particle 
volume  in  capturing  the  initial  behavior  of  the  material  at  higher  Young’s  modulus. 
In  both  cases,  however,  good  correlation  with  the  experimental  data  can  be  seen. 

Throughout  these  comparisons  the  numerical  model  predicts  similar  behavior 
under  plane  strain  compression  as  reported  in  the  experimental  tests.  The  numerical 
model  for  the  196-particle  control  volume  has  not  been  proposed  to  exactly  reproduce 
the  observed  behavior  of  the  plane  strain  compression  of  sand  reported  in  the 
literature.  However,  the  behavior  of  the  algorithmic  specimen  can  only  be  expected 
to  appear  similar  and  consistent  with  the  experimental  results.  The  similarity  of 
the  numerical  and  experimental  responses  demonstrates  the  validity  of  the  numerical 
model  to  predict  plane  strain  responses  of  granular  materials.  Further  extension  of 
the  model  is  now  possible  to  such  areas  as  creep  modeling. 
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Figure  33.  Comparison  between  the  Hertzian 
normal  contact  law  and  the  linearized  contact  law 
for  similar  cylinders. 
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Figure  35.  Deviator  stress  Aa*  versus  axial  strain 
Cl  for  the  49-particle  control  volume  at  different 
Young’s  moduli  and  Osterchelde  sand.  Confining 
stress  =  29.43  N/cm^. 
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—  Experimental  Data 
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Figure  36.  Initial  deviator  stress  Aa*  versus  axial 
strain  ei  for  the  49-particle  control  volume  at 
different  Young’s  moduli  and  Osterchelde  sand. 

Confining  stress  =  29.43  N/cm^. 
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Figure  37.  Deformed  configuration  for  the  49- 
particle  control  volume.  Confining  stress  = 

29.43  N/cm^;  particle  radius  r  =  0.0092  cm. 
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Figure  38.  Deviator  stress  Aa*  versus  axial  strain 
€i  for  the  49-particle  control  volume  at  different 
Young’s  moduli  and  Karlsruhe  sand.  Confining 

stress  =  19.62  N/cm^. 
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Figure  39.  Initial  deviator  stress  Aa*  versus 
axial  strain  ei  for  the  49-particle  control  volume 
at  different  Young’s  moduli  and  Karlsruhe  sand. 

Confining  stress  =  19.62  N/cm^. 


Figure  41.  Initial  configuration  for  the  196- 
particle  control  volume.  Coefficient  of  uniformity 
Cu  —  1.71;  coefficient  of  curvature  Cc  —  0.90; 
mean  particle  radius  =  0.0092  cm. 
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Figure  42.  Sensitivity  analysis:  ramp  function  e. 
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Figure  43.  Sensitivity  analysis:  maximum  normal 
contact  force  Pmax^ 
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Figure  44.  Deviator  stress  Aa*  versus  axial  strain 
e\  for  the  196-particle  control  volume  at  different 
Young’s  moduli  and  Osterchelde  sand.  Confining 

stress  =  29.43  N/cm^. 
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Figure  45.  Initial  deviator  stress  Aa*  versus  axial 
strain  ei  for  the  196-particle  control  volume  at 
different  Young’s  moduli  and  Osterchelde  sand. 

Confining  stress  =  29.43  N/cm". 
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Figure  46.  Deformed  configuration  for  the  196- 
paxticle  control  volume.  Confining  stress  = 

29.43  N/cml 
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Figure  47.  Deviator  stress  Act*  versus  axial  strain 
ei  for  the  196-particle  control  volume  at  different 
Young’s  moduli  and  Karlsruhe  sand.  Confining 

stress  =  19.62  N/cm^. 
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Figure  48.  Initial  deviator  stress  Aa*  versus 
axial  strain  ei  for  the  196-particle  control  volume 
at  different  Young’s  moduli  and  Karlsruhe  sand. 

Confining  stress  =  19.62  N/cm  . 
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Figure  49.  Deformed  configuration  for  the  196- 
particle  control  volume.  Confining  stress  = 

19.62  N/cml 


130 


Chapter  7 

Fluid  Flow  Model 


7.1.  Introduction 

Since  the  dry  model  discussed  in  the  previous  chapters  is  capable  of  predicting 
the  volume  changes  resulting  from  simple  shearing,  we  can  use  it  to  predict  the  pore 
pressure  changes  in  fully  saturated  granular  assemblies  as  a  result  of  changes  in  the 
soil’s  microstructure.  In  reality,  the  pore  pressures  in  a  saturated  granular  assembly 
can  vary  microscopically  in  the  same  random  fashion  as  the  particle-to-particle 
contact  stresses  because  of  local  variations  in  the  volume  of  voids  between  particles. 
For  example,  water  squirted  between  particle  contacts  could  be  one  source  of  such 
local  variation  in  the  pore  water  pressure  [52].  However,  the  overall  pore  pressures 
on  the  macroscale  level  must  be  consistent  with  the  overall  macroscopic  motion.  In 
this  study,  we  assume  that  the  pore  pressure  variable  is  a  macroscopic  quantity  that 
is  homogeneous  throughout  the  cell  of  interest.  A  consequence  of  this  assumption 
is  that  the  presence  of  water  only  serves  to  impose  a  macroscopic  constraint  on  the 
overall  volume  change  of  the  particle  assembly,  but  does  not  directly  impact  the 
local  motion  of  the  individual  particles.  The  following  discussion  elaborates  this 
latter  point. 
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7.2.  Balance  Laws 


Let  &*  be  the  tensor  of  overall  imposed  stresses  on  the  granular  assembly,  & 
be  the  overall  stresses  arising  from  the  particle-to-particle  contacts,  and  0  be  the 
homogeneous  macroscopic  pore  pressure  in  the  assembly;  then  the  eflfective  stress 
equation  reads 

a*  =  &  +  ei  (7.1) 

where  1  is  the  Kronecker  delta.  Now,  let  us  consider  the  influence  of  undrained 
deformation  on  the  overall  deformation  represented  by  the  strain  tensor  6.  Since 
O’  =  ^(e)  and  fluid  flow  is  inhibited,  then  we  can  use  the  constitutive  equation 
0  =  Aty  tr(€)/n,  where  n  is  the  porosity  of  the  assembly,  to  obtain 

a-*  =  d-(e)  +  — tr(e)l  (7.2) 

Tt 

where  A^,  is  the  bulk  modulus  of  the  water  phase.  Now,  if  the  water  phase  itself 
is  assumed  to  be  incompressible,  then  A^,  ^  1,  which  yields  tr(c)  0,  i.e.,  the 
deformation  becomes  macroscopically  volume-preserving.  Note  that  this  constraint 
impacts  the  motion  of  the  particles  collectively  rather  than  individually.  Thus,  the 
volume  constraint  due  to  the  presence  of  water  may  be  viewed  either  as  a  factor  that 
alters  the  overall  stress-strain  behavior  of  the  granular  assembly,  or  as  a  driving  force 
that  causes  a  change  in  the  pore  water  pressure  due  to  the  assembly’s  tendency  to 
change  in  volume. 

An  integral  equation  for  balance  of  momentum  can  be  developed  in  terms  of 
the  overall  stress  tensor  a*.  Balance  of  momentum  over  the  entire  total  domain 
with  volume  U  and  outer  surface  dU  takes  the  form 

f  poGdV+  f  &*  ndA  =  0  (7.3) 

Ju  JdU 

where  po  is  the  reference  mass  density  of  the  infinitesimal  volume  element  dV^  G  is 
the  gravity  acceleration  vector,  and  n  is  the  outward  unit  normal  to  the  surface  dA, 
In  the  context  of  particulate  mechanics,  the  periodic  cell  V  defined  in  the  previous 
chapters  now  takes  the  meaning  of  the  macroscopic  differential  volume  element  dV. 


Assuming  that  both  the  solid  grains  and  fluids  are  incompressible,  balance  of 
mass  for  the  solid-water  mixture  takes  the  form  (see  [53]) 

div  V  -|-  div  5;  =  0  (7.4) 

where  v  is  the  macroscopic  intrinsic  velocity  of  the  solid  phase,  v  is  the  macroscopic 
superficial  Darcy  velocity,  and  div  is  the  spatial  divergence  operator.  For  undrained 
condition,  v  =  0  since  the  solid  and  fluid  phases  move  as  one  body.  Consequently, 
this  condition  gives  rise  to  the  constraint  divr  =  tr(e)  =  0,  which  is  precisely 
recovered  from  (7.2)  when  A^,  ^  1. 

A  mathematical  formulation  casting  (7.3)  and  (7.4)  into  the  framework  of 
a  finite  element  (FE)  continuum  model  is  presented  in  [53].  The  standard  FE 
solution  methodology  entails  dividing  the  domain  into  smaller  finite  elements  and 
evaluating  the  FE  matrix  counterparts  of  equations  (7.3)  and  (7.4).  The  domain 
integrals  are  then  evaluated  numerically  by  Gauss  integration  rule  to  form  a  FE 
consolidation  model.  The  numerical  formulation  of  the  FE  consolidation  model  is 
well  documented  and  is  beyond  the  scope  of  this  report  (see  e.g.  [54]).  However,  the 
link  to  particulate  mechanics  can  be  established  through  the  macroscopic  response 
of  each  Gauss  integration  point,  which  may  be  evaluated  from  the  overall  response 
of  the  particle  assembly  according  to  the  theory  presented  earlier  in  this  report.  In 
the  following  section  we  shall  assume  that  such  link  has  already  been  established, 
and  describe  a  numerical  example  which  illustrates  the  impact  of  fluid  flow  on  a 
continuum  material  model  exhibiting  a  dilatant  behavior. 

7.3.  Numerical  Example 

The  objective  of  this  numerical  example  is  to  study  the  impact  of  fluid 
flow  on  the  deformation  behavior  of  granular  assemblies  treated  as  a  continuum 
material.  Here,  we  consider  a  two-dimensional  plane-strain  FE  mesh  shown  in 
Fig.  50.  The  elements  are  Q9P4  (9-noded  Lagrangian  for  displacements  and  4- 
node  bilinear  for  pore  pressures)  in  which  the  material  behavior  is  given  by  the 
Drucker-Prager  plasticity  model  with  associative  flow  rule.  This  material  model  is 
phenomenological  and  not  derived  from  the  particulate  mechanics  theory.  However, 
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since  the  objective  of  this  example  is  to  study  the  general  behavior  of  the  continuum 
model,  such  substitution  is  not  critical  to  the  results  of  the  formulation.  Ideally,  a 
full  numerical  analysis  would  have  entailed  two  levels  of  modeling,  the  first  being 
the  micromechanical  level  and  the  second  being  the  continuum  level;  however,  the 
computations  engendered  by  this  type  of  analysis  are  extensive  and  prohibitively 
expensive,  and  so  the  first  level  of  modeling  will  be  suppressed  altogether  and  will 
be  replaced  with  a  phenomenological  constitutive  model. 

Figure  50  also  shows  all  the  boundary  conditions  and  imposed  deformation. 
The  elastic  material  parameters  are:  Young’s  modulus  E  =  5  x  10^  kPa,  and 
Poisson’s  ratio  v  =  0.1;  regular  elements  have  cohesion  c  =  80  kPa,  friction  angle 
^  =  35®,  and  hardening  parameter  H'  =  10  kPa;  weak  elements  are  placed  on  the 
side  of  the  mesh  of  Fig.  50  to  trigger  localized  deformation,  and  are  defined  by  a 
material  having  cohesion  of  c  =  30  kPa,  friction  angle  of  ^  =  15®,  and  hardening 
parameter  of  H'  =  0.  All  boundary  pore  water  pressures  are  assumed  zero. 

Figure  51  shows  plots  of  average  normal  stress  in  the  X2-direction  (d’22)  versus 
the  imposed  average  nominal  strain  (€22),  assuming  the  response  to  be  fully  drained 
(no  excess  pore  water  pressure).  The  deformation  is  applied  in  increments  of 
0.0625  m,  or  1.25%,  with  a  final  compressed  deformation  of  20%.  The  mesh  with 
weak  elements  substituted  at  the  side  begins  to  yield  and  reaches  a  plateau  at  stresses 
lower  than  those  attained  without  weak  elements;  hence,  the  stress-strain  curve  for 
the  mesh  with  weak  elements  plots  on  the  lower  side  relative  to  the  curve  for  the 
mesh  without  weak  elements.  This  phenomenon  is  a  standard  result. 

Figure  52  shows  plots  of  average  normal  stress  in  the  X2-direction 
versus  the  imposed  average  nominal  strain  (€22)  for  meshes  without  weak  elements, 
assuming  the  responses  to  be  fully  drained  and  fully  undrained.  For  the  fully 
undrained  case,  consolidation  is  turned  on  at  the  end  of  the  last  loading  phase, 
with  the  total  imposed  deformation  held  fixed  during  this  transient  period  of  pore 
pressure  dissipation.  Note  that  the  stresses  produced  in  the  mesh  for  the  undrained 
case  are  greater  than  those  reached  for  the  drained  case  because  the  pore  water 
pressure  is  negative  due  to  a  dilating  mesh.  Consequently,  the  inclusion  of  the  fluid 
flow  effect  “strengthens”  the  mesh  for  this  case  dilating  material.  Again,  this  is  a 


standard  result. 

While  the  results  presented  in  this  chapter  are  fairly  standard,  it  sheds  hght 
onto  the  possible  connection  between  the  micromechanical  theory  presented  in  this 
report  and  the  e^dsting  continuum  models  available  for  granular  materials.  In 
principle,  this  connection  should  be  clear  from  the  outset  as  elaborated  in  this 
numerical  example.  The  modeling  approach  presented  in  this  report  is  not  the 
same  as  the  many  particulate  mechanics  models  available  in  the  literature,  in  which 
the  micro-  and  the  macromechanical  aspects  of  the  problem  are  combined  directly 
to  form  a  single  particle  assembly  that  is  subjected  to  a  non-homogeneous  mode 
of  deformation.  The  latter  approach  may  have  some  difficulty  tracing  the  actual 
effect  of  material  behavior  and  that  of  the  nonhomogeneous  boundary  condition. 
Computationally,  the  modeling  approach  presented  in  this  study  is  amenable  to 
parallel  computation,  and  with  the  current  advances  in  computational  hardware 
and  software  modeling  tools,  the  prognosis  for  a  full  analysis  entailing  two  levels  of 
modeling  is  indeed  positive. 


Figure  51.  Resulting  average  stress  ^22  versus  imposed  average  strain  €22  demonstrating 
lower  stress  plateau  reached  when  two  weak  elements  are  incorporated  at  side  of  mesh. 


Figure  52.  Resulting  average  stress  (732  versus  imposed  average  strain  €22  demonstrating 
“strengthening”  effect  due  to  negative  pore  pressures  generated  in  a  dilating  material. 


Chapter  8 

Summary  and  Conclusions 


A  methodology  has  been  presented  for  determining  the  overall  response  of 
a  granular  material  based  on  the  fundamental  mechanism  of  particle-to-particle 
interaction.  The  mathematical  foundation  rests  on  a  particulate  mechanics 
description  of  the  behavior  of  two  contacting  particles.  With  this  micromechanical 
framework,  a  numerical  algorithm  solves  the  problem  of  a  prescribed  deformation 
history  for  an  assembly  of  particles  representing  a  point  in  a  granular  material. 
To  make  the  overall  assembly  response  independent  of  the  imposed  boundary 
displacements,  the  formulation  employs  the  notion  of  a  repeating  cell.  An  important 
feature  of  the  mathematical  formulation  lies  in  its  generation  of  the  micromechanical 
responses  quasi-statically,  and  not  dynamically.  These  fundamental  features  make 
the  model  amenable  to  incorporation  into  conventional  finite  element  codes. 

To  assess  the  performance  of  the  model,  two-dimensional  plane-strain 
simulations  with  regular  and  random  initial  packing  of  circular  disks  were  performed. 
The  model  captures  some  of  the  most  important  features  of  granular  material 
behavior  such  as  anisotropy,  hardening,  and  softening  responses.  In  general,  the 
numerical  stability  of  the  solution  depends  upon  the  character  of  the  imposed 
overall  motion — convergence  of  the  iteration  is  easy  whenever  the  shearing  is 
preceded  by  isotropic  compression,  but  is  difficult  whenever  the  imposed  motion 
lacks  the  isotropic  compression  needed  to  keep  the  microstructure  from  collapsing. 


Interestingly,  this  characteristic  of  the  numerical  model  reflects  the  property  of  the 
subject  prototype — that  the  stress  response  of  granular  media  depends  strongly 
upon  the  imposed  confining  load. 

In  a  natural  extension  of  the  deformation  driven  model,  a  methodology  has 
been  developed  for  determining  the  overall  deformation  response  of  a  granular 
material  to  an  imposed  stress  history.  This  formulation  allows  finite  motions  of 
the  individual  particles,  but  the  macroscopic  material  response  is  postulated  based 
on  infinitesimal  theory.  An  important  by-product  of  the  mathematical  formulation 
is  the  overall  tensor  of  material  moduli  derived  from  the  consistent  linearization 
of  the  elasto-plastic  constitutive  equation.  Owing  to  the  complexity  associated 
with  the  implementation  of  the  exactly  linearized  algorithm,  an  alternative  secant 
approximation  has  been  implemented  that  retains  the  essential  properties  of  the 
exact  description.  These  moduli  are  then  used  to  predict  the  onset  of  localized 
deformation  in  granular  materials  on  the  macro-scale  level.  Again,  a  unique  feature 
of  the  model  is  its  use  of  first  principles  to  predict  the  stability  behavior  of  an 
assembly  of  discrete  particles. 

Numerical  simulations  with  uniform  Jis  well  as  non-uniform  assemblies  of  two- 
dimensional  circular  disks  demonstrate  the  model’s  capability  to  capture  some 
of  the  most  important  features  of  granular  material  behavior.  The  stress-driven 
format  of  the  numerical  algorithm  allows  the  model  to  capture  naturally  such 
mechanical  responses  as  structural  anisotropy,  pressure  dependency,  and  volume 
change  characteristics  of  granular  materials.  Thus,  a  fundamental  understanding 
of  the  behavior  of  particulate  materials  is  available  based  on  the  interaction  of  the 
constituent  particles. 


Appendix  A.  Virtual  contact  separation 


To  calculate  A°'  in  equation  (3.8)  consider  two  rigid  control  volume  particles 
A  and  B  in  contact  at  c.  Let  the  virtual  displacement  field  produce  a  relative 
displacement  at  contact  c  given  by 

A"  =  ,  (Al) 

where  represents  the  virtual  displacement  of  c  from  the  perspective  of  particle  A. 

Since  the  particles  are  assumed  rigid,  the  relative  motion  of  the  particle  contact 
can  be  written  in  terms  of  a  translation  and  a  rotation.  For  example,  particle  A’s 
contact  motion  can  be  written  as 

•  (*"  -  x^)  ( A2) 

where  =  —w^'^  represents  the  rotation  of  contact  c  relative  to  A,  and 
denotes  the  displacement  of  the  centroid  of  A. 

Substituting  equation  (A.2)  and  a  similar  equation  for  particle  B  into  equation 
(A.l)  gives 

A^  =  u^  -u^  +  .  {x^  -  x^)  -  -  x^) .  (A3) 


The  displacement  of  the  particle  centroids  can  be  defined  from  the  series 
expansion  about  the  common  contact  point  c.  Assuming  that  the  displacements 
and  and  the  rotations  and  conform  to  some  smooth  field  u  and  w 
and  that  =  u{x^)  represents  the  value  of  u  at  the  contact  c,  the  contact  rotations 
become  —  w{x^) .  The  centroid  displacements  can  be  written  as 


dui{x*^) 

dxk 


i^k  -  ^k)  + 


(A4) 


and 


Uj  —  U  +  (ijj.  Xj.)  +  .  .  .  . 


Subtracting  (A.4)  from  (A.5)  and  substituting  into  (A.3)  gives 


dui{x‘^) 


A  n 


(A5) 


If  the  tensor  <l>  is  identified  as 


<f>ik  = 


dui{x*^) 

dxk 


-  Wikix*") , 


then  to  a  first  order  approximation  the  virtual  separation  reads 


which  is  identical  to  the  field  equation  (3.10)  specified  for  the  contact 


Appendix  B.  Localization  function 


For  a  path-dependent  material  response,  the  constitutive  equations  are 
integrated  incrementally.  When  the  deformation  reaches  the  point  where  the 
localization  function 

/(n)  =  det[A(n)]  (B.l) 

becomes  negative  or  zero,  localization  is  said  to  have  begun.  To  predict  this  onset 
of  localization  the  minima  of  the  function  /  must  be  computed  at  the  end  of  every 
deformation  increment. 


In  the  two-dimensional  case,  the  acoustic  tensor  A(n)  takes  the  following  form 


A{n) 


n\Kun  +  nin2K\n2 

■¥n2n\K\2n  +n\K\2n 

wfA'mi  +  nin2K\2\2 
.-\-n2n\K22\\  ■^n\K22l2 


n\Kn\2 n\n2Kii22 

■\-n2n\K\2\2  +  'n\K\222 

n\K\2\2  +  n\n2K\222 
+n2n\K22\2  +  ”2^2222 . 


(S.2) 


where  ni  =  cosO  and  n2  =  sin0,  and  the  angle  6  defines  the  orientation  of  the 
discontinuity  planes.  The  localization  function  yields 


/(n)  =  aonf  -|-  ainin2  -|-  a2n\n\  -f  aanin^  H-  04112 


(B.3) 


where 

Oo  =  -^^1111-^^1212  —  -h^lll2^1211 

Ol  =  -h'llllJ'i'1222  +  ■K’llll-ft’2212  —  -h"lll2f^2211  —  -^"1122^^1211 

02  =  -K’iiii.K'2222  +  -^1112^1222  +  -^^121 1-^^2212 

—  -fiLll22f^l212  —  ^1122^2211  “  fii^l212^2211 
03  =  K1112K2222  +  Ki2nK2222  -  Kn22K22\2  —  Ki222K22n 
04  =  K1212K2222  —  K1222K22U 

With  the  change  of  variables  x  =  tan  0,  equation  (B.3)  can  be  written  as 

f{x)  =  040:^  -t-  a^x^  -f  a2X^  -f-  oix  +  oq. 


{BA) 


(B.b) 


The  minima  of  the  function  f{x)  occur  at  the  roots  of  the  cubic  polynomial  f'{x) 
which  can  be  analytically  found  using  Carden’s  formulae.  Positive  minima  indicate 
that  localization  has  not  yet  developed  while  a  negative  or  zero  minimum  signals 
the  onset  of  localization. 


Appendix  C.  Strain-driven  problem:  tangent  operator  derivatives 

Unless  otherwise  specified,  all  quantities  in  this  section  have  implied  subscripts 
of  “n  +  1”  for  time  step  and  superscripts  of  ‘i’  for  iteration  step. 

Taking  particle  A’s  point  of  view,  we  have  the  following  preliminary  results: 

A  :=  n'(d^)  =  -n'(d^)  =  -(I  -  nn*)/||  I  i| ;  ((7.1) 

n'(0^)=n'(0^)  =  O,  ((7.2) 


where  n  =  n^.  Note  that  A  =  A^,  or  Aij  =  Aji.  Other  necessary  derivatives  are  as 
follows.  For  the  normal  indentation,  one  has 


S^{d^)  =  -6'{d^)  =  n-, 

(C.3) 

S'{e^)  =  S'{0^)  =  O. 

((7.4) 

One  can  verify  that  for  the  incremental  slip. 

A7'(d^)  =  -A7'(<i^)  =  (r^  +  r^)A^‘^'(d'^)  i 

((7.5) 

A7'(^"‘)  =  -  ; 

((7.6) 

Ay  (0^)  =  -  , 

((7.7) 

where 

A0^'(d^)  =  (sec  A^^n„  x  63  +  tanA^^n)/||/  || ; 

((7.8) 

A0^  =  sin“^(e3  •  n„  x  n) . 

((7.9) 

Also,  the  derivative  of  the  ramp  function  becomes 

(  If  ^i»+i  ^ 

HUS)  =  {l/e,  ifO<^„+i<e; 

1  0,  if  <  0. 

((7.10) 

Now,  assuming  that  particles  A  and  B  are  in  contact,  one  finds  that 


.  fW'")  =  -A(d^)  =  +  H.kf,]  ; 

/«(«•’)  =  A(»®)=0. 


(ail) 

((7.12) 


If  the  contact  mode  is  ‘stick,’  the  corresponding  derivatives  of  the  tangential  force 
become 

/r(<l'')  =  =  H',kT^i  g(d*)  +  H,hr  Ai\d*) ,  (C.13) 

/!■(«•’)  =  (C.14) 

and  if  the  contact  mode  is  ‘slip,’  the  derivatives  of  the  tangential  force  are 

/f  (d'*)  ^  -/f  (d«) ,  (C.15) 


/'.(d-’j  = 


kr 


H.kr  +  H' 


H'HeA'y'id^)  +  fe/r'(A7  -  A7P) 

-  {•)H,kNtan<l>{HX+i+H.)]  S'{d^) 


.(.A.  KH’kT  A 

where  (•)  =  (/j’)n  +  HekrAy. 

Constructing  the  local  force  gradients,  one  obtains 

=  I 


\{2x2)  ’ 


(2x1) 

Since  f  =  R-  the  chain  rule  yields 

f'(d^)  =  -fid^)  =  R .  J^{d^)  +  R!{d^)  •  :F. 


,  ((7.16) 
((7.17) 


(C.18) 

((7.19) 

((7.20) 


The  last  term  on  the  right-hand  side  of  ((7.20)  has  the  explicit  form 


R'{d^) .  J='  = 


(An  In  -  ^12/r) 
(Au/n  +  ^nfr) 


{AnfN  -  M2fT) 
{M2fN  +  ■di2/r) 


((7.21) 


in  which  the  ylij’s  are  given  in  ((7.1).  Finally, 

f'{e^)  =  (^)  f{e^)  =  R .  -i-  R!{e^)r  =  Me^)  n  x  ca .  ((7.22) 

This  completes  the  derivatives  necessary  to  construct  the  gradient  f^'{d^). 


Appendix  D.  Stress-driven  problem:  tangent  operator  derivatives 


Consider  contact  element  ‘e’  to  represent  the  contact  between  particle  A  and  B. 
Unless  otherwise  specified,  all  quantities  in  this  section  have  implied  superscripts  of 
‘A:’  for  iteration  step.  Taking  particle  A’s  point  of  view  to  write  the  force  derivatives 
in  (5.46)  and  (5.47)  can  be  written  as 


dd^  L  dd^  J  (3x3) 

(B.l) 

dd^  ~  L  ddA  de^  j(3^3) 

(D.2) 

• 

where 

(D.3) 

• 

and 

1  ^  J(3xl)  I*'  J(3xl) 

(DA) 

The  following  preliminary  results  extend  from  Appendix  C: 

• 

4  drift  II  drift  ^  drift 

■“  dd''  ~  ^  ~ 

=  0  ;  (Z>.5) 

• 

where  and  Un+i  =  •  Note  that  An  =  A*-  and  An+i  =  A\^-y . 

Other  necessary  derivatives  are  as  follows.  For  the  normal  indentation,  one  has 


dSn+l  dSn+l  dSn+1  „ 

adA  =  adA  “  ’  d0A  -  ®- 

(D.l) 

m 

For  the  incremental  slip,  one  can  write 

dA7  dA7ft  dAya 
dd^  “  dd^  dd^  ’ 

(D.8) 

• 

and 

dAy  _  dA7ft  dAya 
dd^  ”  dd^  dd^  ’ 

(D.9) 

145 


where 


and 


dd^ 

dA'fn 

dd^ 

dA'fa 

ddA 


=  { 


4-  r 


dAO^ 

f  d~dA  ’ 


A 


=  secA0^  n„  x  63  +  tan  ng  /||Zfi|| , 


=  0, 


rA  ^b\ 


dA0^ 


(D.IO) 

(D.ll) 

(n.i2) 

(£>.13) 


=  ( 

=  +  r^)  [sec  A0?  €3  x  n„+i  +  tan  A^f  /||/ft|| 

+  +  r'®)  [sec  A0f  ng  x  63  +  tan  A^f  nn+i  /l|Z„+i||  ,(£>.14) 

_  [sec  A^?  fin  X  €3  +  tan  A^?  n„+i]  /||ln+i||  •  (£>.16) 

Similarly,  considering  the  rotation  effects  gives 


dA'/A 

ddA 


dA'f  dAjn  dAjA 
dOA  ~  dQA  dQA  ’ 

where 


dA-^n 

d$A 

dA'in 

doA 


The  derivative  of  the  ramp  function  becomes 


{0,  if  5„+i  >  c ; 

1/e,  if0<<5„+i<c; 

0,  if  <  0 . 


Now,  the  derivatives  of  the  normal  force  can  be  written  a.s 


dfN  dfN 


ddA 

dfN 

aeA 


[H^k^Sn+i  +  HekN]  n„+i , 


(£>.17) 

(Z?.18) 

(£>.19) 

(£>.20) 

(£>.21) 

(£>.22) 


If  the  contact  mode  is  ‘stick,’  the  corresponding  derivatives  of  the  tangential  force 
become 


Bd^  '  Bd^'  Bd* 

(D.23) 

dd^  dd^  dd^ 

(D.24) 

^IL  —  —r^Hekr 

QQA  ^ 

(Z).25) 

If  the  contact  mode  is  ‘slip,’  the  derivatives  of  the  tangential  force  are 


1 


dfr  _ _ 

dd^  Hckx  +  H' 


-  A7P) 


—  {•)H(kTkN  tan  <i>{H[8n^\  +  H() 


dd^ 


+  H'H.kr^^ 


dfr  _  1 

dd^  ~  H,kT  +  H'[i 


dd^ 

H[kTH'{Aj  -  A7P) 


—  {•)H(kTk;f  tan  <j>(H'^6n+i  +  Hf)^ 

dA-y 


dSn+i 

dd^ 


+  H'H.kj 


dd^\ ’ 


where  (•)  =  +  HfkxA'j. 

Since  f  =  R  -  the  chain  rule  yields 


Sf  I  SRjd'') 


ad-f  dd^  Bd^ 

Bd^  Bd^  BdJ' 


where 


dr  _ 

9fN 

dfr' 

dd^  ~ 

[dd^  ’ 

dd^. 

dr 

\dfN 

dfr 

ddA  ~ 

[dd^  ’ 

ddA 

T  t 


(2x2) 


(D.26) 


(D.27) 


(I>.28) 

(D.29) 

(D.30) 

(I>.31) 
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The  last  term  on  the  right-hand  side  of  (D.28)  and  {D.29)  has  the  explicit  form 


dR{dA)  _  dR{dA)  _  {Au/n  —  Aufr)  {Au/n  -  ■^22/r) 

dd^  dd^  L  (^12/jv  +  Anfr)  (A22fN  +  Anfr) 


{D.S2) 


in  which  the  Aij's  are  given  in  (/).6).  Finally, 

df  „  dT  .  dR{d^)  ^  dT  _ 

de^  d0^  d9A  ^  ' 


(D.33) 


This  completes  the  derivatives  necessary  to  construct  the  algorithmic  tangent 
operator. 
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Computer  Code  for  Micromechanical  and  Macromechanical  Model 


PROGRAM  MACRO 


....  MAIN  PROGRAM 

IMPLICIT  REAL*8CA-H,0-Z) 

. . . .  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 
DIMENSION  IDIAGC5000) 

DIMENSION  AUPPER(200000) ,ALOWER(200000) ,8(3000) 
DIMENSION  XAC2),XB(2),XL0(2) 

COMMON  /MACRO  /  FSIGC3),SIGC3),EP(3) 

COMMON  /PARTCL/  X(3,1000),NPART,RMIN 
COMMON  /INFO  /  IEN(2, 20000), ID(3, 1000), LMC6, 20000), 
X  NUMEL,NEQ,NSTIFF 

COMMON  /DSDATA/  D(3, 1000), 00(3, 1000), DINCC3, 1000), 

X  DDINC(3,1000) 

COMMON  /FDATA  /  FTN(20000),ALPHA(20000),ICON(20000) 
COMMON  /ELDATA/  SKN,SKT,A,TANPHI,HP 
DATA  ZERO,ONE/0.0D0,1.0D0/ 

0PEN(UNIT=12 , FI LE= • outl ’ , STATUS= ’ UNKNOWN ' ) 
0PEN(UNIT=24, FILE= ' out2 ’ , STATUS= ' UNKNOWN ' ) 
0PEN(UNIT=36, FILE=’ out3 ' , STATUS=’ UNKNOWN ' ) 
OPEN(UNIT=48 , FILE=’ out4' , STATUS=’ UNKNOWN ' ) 
OPEN(UNIT=60 , FILE=’ outs ' , STATUS=' UNKNOWN ' ) 
0PEN(UNIT=51 , FI LE= • input' , STATUS= ’ UNKNOWN ’ ) 

....  INPUT  CONTROL  DATA. 

THE  INPUT  FILE  'input'  TAKES  THE  FOLLOWING  FORM: 

NPART,XCELL,YCELL,SKN,SKT,A,PHI,HP 

X(1),Y(1),RADIUS(1) 


X(NPART) , Y(NPART) , RADIUS(NPART) 

WHERE: 

NPART=NUMBER  OF  PARHCLES  IN  THE  CONTROL  VOLUME 
XCELL=X  DIMENSION  OF  THE  CONTROL  VOLUME. 

YCELL=Y  DIMENSION  OF  THE  CONTROL  VOLUME. 
SKN=NORMAL  SPRING  STIFFNESS. 

SKT=TANGENTIAL  SPRING  STIFFNESS. 

A=PARTICLE  COHESION. 


C  PHI=PARTICLE  FRICTION  ANGLE. 

C  HP=HARDENIN6  PARAMETER. 

C  TANPHI=TAN(PHI). 

C  RADIUS(N)=RADIUS  OF  PARHCLE  N. 

C  X(N)=X  COORDINATE  OF  THE  PARTICLE  N’S  CENTROID. 

C  Y(N)=Y  COORDINATE  OF  THE  PARTICLE  N'S  CENTROID. 

C 

C  THIS  DATA  IS  READ  INTO  THE  ARRAY  X  WHERE  EACH  ELEMENT  HAS 
C  THE  FOLLOWING  MEANING: 

C 

C  X(1,N)  =  X(N) 

C  XC2,N)  =  Y(N) 

C  X(3,N)  =  RADIUSCN). 

C 

C  RMIN  =  MINIMUM  PARTICLE  RADIUS  IN  THE  CONTROL  VOLUME. 

C 

READ(51 , ♦)  NPART.XCELL, YCELL , SKN , SKT, A , PHI , HP 
WRITE(12,999)NPART,XCELL,YCELL,SKN,SKT,A,PHI,HP 
WRITEC6 , 999)NPART,XCELL , YCELL , SKN , SKT, A , PHI , HP 
WRITE(36, 999)NPART, XCELL , YCELL , SKN , SKT, A , PHI , HP 
c  WRITEC12,996) 

RMIN  =  1.0D5 
RMAX  =  0.0D0 
DO  15  1=1,  NPART 

READC51,*)  X(1,I),X(2,I),X(3,I) 

RD  =  X(3,I) 

RMIN  =  DMIN1(RMIN,RD) 

RMAX  =  DMAX1(RMAX,RD) 
c  WRITE(6 , 998)X( 1 , I) , X(2 , I) , XC3 , I) 

15  CONTINUE 

PI=4.0D0*DATANC1.0D0) 

TANPHI=DTANCPHI*PI/180.0D0) 

C 

C....  BOUNDARY  CONDITION  DATA:  FIX  ORIGIN  PARHCLE 
C 

CALL  ICLEAR(ID,3*NPART) 

DO  1  1=1,  3 
ID(I,1)  =  1 
1  CONTINUE 

c  do  2  1=1,  npart 
c  id(3,i)=l 

c  2  continue 
C 

C....  ESTABLISH  EQUATION  NUMBERS 
C 

NEQ  =  0 

DO  10  N=l,  NPART 
DO  10  1=1,  3 

IFCID(I,N).EQ.0)  GOTO  5 
ID(I,N)  =  0 
GO  TO  10 

5  NEQ  =  NEQ  +  1 


ID(I,N)  =  NEQ 
10  CONHNUE 
c  WRITEC12,994) 

c  DO  11  1=1,  NPART 

c  WRITE(12,993)  I,  (10(1, I),  J=l,3) 

c  11  CONTINUE 
C 

C....  INPUT  CONTACT  ELEMENT  DATA 
C  ASSUME  THAT  A  PARTICLE  MAY  INTERACT  ONLY  WITH 
C  PARTICLES  WITHIN  A  RADIUS  EQUAL  TO  4*(LARGEST 

C  PARHCLE  RADIUS). 

C 

XA(1)  =  XCELL 
XA(2)  =  ZERO 
XB(1)  =  ZERO 
XB(2)  =  YCELL 

DET  =  XA(1)*XB(2)-XB(1)*XA(2) 

NUMEL  =  0 

DO  20  NA=1,NPART-1 
DO  60  NB=NA+1, NPART 
C 

RA  =  X(3,NA) 

RB  =  X(3,NB) 

RPLUSR  =  RA  +  RB 
C 

XL0(1)  =  X(1,NB)  -  X(1,NA) 

XL0(2)  =  X(2,NB)  -  XC2,NA) 

AA  =  (XB(1)*XL0C2)-XBC2)*XL0C1))/DET 
BB  =  (XA(2)*XL0(1)-XAC1)*XL0(2))/DET 
SIGNA  =  ZERO 

IF(AA.NE.ZERO)  SIGNA=AA/DABSCAA) 

SIGNB  =  ZERO 

IF(BB.NE.ZERO)  SIGNB=BB/DABSCBB) 

C 

DO  65  1=0,1 
DO  70  J=0,1 
DI  =  DREAL(I) 

DJ  =  DREAL(J) 

XL0(1)  =  X(1,NB)  +  DI*SIGNA*XA(1)  +  DJ*SIGNB*XBC1) 
XL0C2)  =  X(2,NB)  +  DI*SIGNA*XA(2)  +  DJ*SIGNB*XBC2) 
XL0C1)  =  XL0C1)  -  XC1,NA) 

XL0(2)  =  XL0C2)  -  X(2,NA) 

DELTA  =  DSQRT(DOT(XL0,XL0,2)) 
IF(DELTA.LT.4.0D0*RMAX)  GOTO  55 
70  CONHNUE 
65  CONTINUE 
GOTO  60 
55  CONHNUE 

NUMEL  =  NUMEL  +  1 
IEN(1, NUMEL)  =  NA 
IEN(2, NUMEL)  =  NB 
60  CONHNUE 


CONHNUE 


20 
C 

C....  INITIALIZE  THE  HARDENING  ARRAY  TO  PARTICLE  COHESION  AND 
C  THE  CONTACT  ARRAY  TO  ALL  PARTICLES  IN  CONTACT  WITH  EACH 
C  OTHER  (ie  ICON(NUMEL)=l). 

C 

DO  25  1=1,  NUMEL 
ALPHACI)=A 
ICONCI)=l 
25  CONTINUE 
C 

C....  LOCALIZE  ID  ARRAY 
C 

DO  40  K=l,  NUMEL 
DO  30  J=l,  2 
NN  =  lENCJ.K) 

DO  30  1=1,3 

LM(3*(J-1)+I,K)  =  IDCI,NN) 

30  CONTINUE 
40  CONTINUE 
C 

C....  COMPUTE  DIAGONAL  ADDRESSES 
C 

NSHFF  =  1 
IDIAG(l)  =  1 
IF(NEQ.EQ.l)  GOTO  100 
DO  50  1=2,  NEQ 

IDIAGCI)  =  IDIAGCI-1)  +  I 
50  CONTINUE 

NSTIFF  =  IDIAG(NEQ) 

C 

C....  CALL  DRIVER  PROGRAM 
C 

100  CALL  MACROI(XCELL,YCELL,AUPPER,ALOWER,B,IDIAG) 

C 

C....  PRINT  RESULTS 
C 

993  F0RMAT(2X,I3,9X,3I4) 

994  FORMATC/, ’ID  ARRAY: ',/,’PARnCLE',3X, 'GLOBAL  EQ.  NUMBER') 
996  FORMATC// , 10X , ' X ' , 15X , ' Y ' , 14X , ' RADI US ' ) 

998  F0RMAT(F16.8,F16.8,F16.8) 

999  FORMATC 'NUMBER  OF  PARHCLES  =  ',13,/, 

X  'WIDTH  OF  THE  CONTROL  VOLUME:  XCELL  =  ',F16.8,/, 

X  'HIEGHT  OF  THE  CONTROL  VOLUME:  YCELL  =  ',F16.8,/, 

X  'NORMAL  SPRING  SHFFNESS  =  ',E20.8,/, 

X  'TANGENTIAL  SPRING  STIFFNESS  =  ',E20.8,/, 

X  'PARTICLE  COHESION  =  ',E20.8,/, 

X  'PACnCLE  FRICnON  ANGLE  =  ',F16.8,/, 

X  'HARDENING  PARAMETER*  ',E20.8,/) 

C 

CLOSECUNIT=12) 

CLOSECUNIT=24) 


CL0SECUNIT=36) 

CL0SECUNIT=48) 

CLOSECUNIT=60) 

CL0SE(UNIT=51) 


SUBROUTINE  MACROI(XCELL,YCELL,AUPPER,ALOWER,B,IDIAG) 

C 

C....  INVERSE  PROBLEM  DRIVER  PROGRAM 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

C 

C....  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 
C 

DIMENSION  XAOLD(2) , XBOLD(2) , EPSI(3) , EPSC3 , 1000) , DELSIG(3, 1000) 
DIMENSION  EK(6 , 6) , P(6) , XA(2) , XB(2) , DSIG(3) , SK(3 , 3) 

DIMENSION  AUPPER(l) .ALOWERCl) ,B(1) , IDIAGCl) 

DIMENSION  UNITY(3 , 3) , SIG1(3) , SIG2(3) , SIGTC3) 

COMMON  /MACRO  /  FSIG(3),SIG(3),EP(3) 

COMMON  /PARTCL/  X(3,1000),NPART,RMIN 

COMMON  /INFO  /  IENC2, 20000), ID(3, 1000), LMC6, 20000), 

X  NUMEL,NEQ,NSTIFF 

COMMON  /DSDATA/  D(3, 1000), DD(3, 1000), DINC(3, 1000), 

X  DDINCC3,1000) 

COMMON  /FDATA  /  FTN(20000),ALPHAC20000),ICONC20000) 

COMMON  /ELDATA/  SKN,SKT,A,TANPHI,HP 
DATA  ZERO,ONE/0.0D0,1.0D0/ 

C 

CALL  CLEAR(UNITY,9) 

DO  1  1=1,  3 

UNITY(I,I)=ONE 
1  CONTINUE 
C 

C....  CLEAR  MACRO  QUANTITIES:  STRESS  AND  STRAIN 
C 

CALL  CLEAR(FSIG,3) 

CALL  CLEAR(SIG,3) 

CALL  CLEAR(EP,3) 

CALL  CLEAR(SIG1,3) 

CALL  CLEAR(SIG2,3) 

C 

C....  CLEAR  MICRO  QUANTITIES:  TANGENTIAL  FORCE,  DISPLACEMENT, 

C  AND  DISPLACEMENT  INCREMENT. 

C 

CALL  CLEAR(FTN,NUMEL) 

CALL  CLEAR(D,3*NPART) 

CALL  CLEAR(DINC,3*NPART) 

CALL  ICLEARCICON.NUMEL) 

C 

C....  INIHAL  CONTROL  VOLUME  CONFIGURATION 
C 

XA(1)  =  XCELL 
XA(2)  =  ZERO 
XBCD  =  ZERO 
XB(2)  =  YCELL 
C 

C....  ANALYSIS  OPTIONS 


ICODE=0  =>  STRAIN  CONTROL 
IC0DE=1  =>  STRESS  CONTROL 

IC0DE=2  =>  2  TIMESTEPS  STRAIN  CONTROL;  STRESS  CONTROL 
I CODE  =  2 

IF(ICODE.EQ.l)  GOTO  7000 


STRAIN  CONTROL 


. . . .  INPUT  MACROSCOPIC  STRAIN  INCREMENTS  (EPS) 

A(k-l)=eps(I,l)  A(k)=eps(I,2) 
WRITEC6,980) 

WRITE(36,980) 

NSTEP  =  2 
f  ntnax=0 . 0d0 

CALL  CLEAR(EPS,3*NSTEP) 

DO  10  1=1,  NSTEP 

EPSCl.I)  =  -8.500d-3 
EPSC2,I)  =  -8.500d-3 
EPSC3,I)  =  1.00d-5 

IF(I.GE.2)  THEN 

EPS(1,I)  =  -8.600d-3 
EPS(2,I)  =  -8.600d-3 
EPS(3,I)  =  2.00d-5 
ENDIF 
10  CONTINUE 

. . . .  TIME  STEP  LOOP 

DO  1000  NS=1, NSTEP 
WRITEC6,982)  NS 
WRITEC12,982)  NS 
WRITEC24,982)  NS 
WRITE(36,982)  NS 

....  UPDATE  CONTROL  VOLUME  CONFIGURATION 
DO  15  1=1,3 

EP(I)  =  EP(I)  +  EPS(I,NS) 

EPSI(I)  =  EPS(I,NS) 

15  CONTINUE 

DO  150  1=1,2 
XAOLD(I)=XA(I) 

XBOLD(I)=XBCI) 

150  CONTINUE 


C....  CALL  STRAIN  DRIVER  PROGRAM 
C 

CALL  MI CROC EPSI , XA , XB , XAOLD , XBOLD , NS , NSTEP , DET, 

X  AUPPER,ALOWER,B,IDIAG,fnniax) 

C 

C....  UPDATE  TANGENHAL  CONTACT  FORCES/CALCULATE  MACRO  STRESS 
C 

CALL  CLEAR(SIG,3) 

DO  2000  N=1,NUMEL 
NA  =  lENCl.N) 

NB  =  IEN(2,N) 

ITASK  =  1 

IF(ICODE.EQ.2)  ITASK=2 

CALL  CNTACTCN , NA , NB , P , EK , XA , XB , DET , ITOUCH , ITASK , XAOLD , 

X  XBOLD, fnmax.npl as) 

2000  CONTINUE 
C 

XAxXB  =  DABS(XA(1)*XB(2)  -  XA(2)*XB(1)) 

AREA  =  XAxXB 
SIG(l)  =  SIG(1)/AREA 
SIG(2)  =  SIG(2)/AREA 
SIG(3)  =  SI6C3)/AREA 
C 

IF(ICODE.NE.2)  GOTO  4999 

XA(1)=XCELL 

XA(2)=ZERO 

XB(l)=ZERO 

XB(2)=YCELL 

IF(NS.EQ.l)  THEN 

CALL  MOVECSIGl,SIG,3) 

GOTO  5001 
ELSE 

CALL  MOVE(SIG2,SIG,3) 

GOTO  5001 
ENDIF 
C 

C....  UPDATE  COORDINATES  OF  PARTICLE  CENTROIDS 
C 

4999  CONHNUE 

DO  5000  1=1, NPART 
X(1,I)  =  X(1,I)  +  DINC(1,I) 

X(2,I)  =  X(2,I)  +  DINC(2,I) 

C 

C....  CHECK  FOR  IMAGE  PARHCLE 
C 

AA  =  (XB(2)*X(1,I)-XB(1)*XC2,I))/DET 
BB  =  (XA(1)*X(2,I)-XA(2)*X(1,I))/DET 
IF(AA.GE.ZERO.AND.AA.LE.ONE.AND.BB.GE.ZERO.AND.BB.LE.ONE) 
X  GOTO  5000 
C 

C....  CONVERT  IMAGE  PARTICLE  TO  REAL 


IF(AA.LT.ZERO)  AA=AA+ONE 
IF(AA.6T.0NE)  AA=AA-0NE 
IF(BB.LT.ZERO)  BB=BB+ONE 
IF(BB.GT.ONE)  BB=BB-ONE 
XCl.I)  =  AA*XA(1)  +  BB*XB(1) 

XCZ.I)  =  AA*XAC2)  +  BB*XB(2) 

5000  CONTINUE 

....  OUTPUT 

5001  CONTINUE 
PRINCIPAL  STRESSES 

TEMP=DSQRT((SIGC1)-SIG(2))*(SIG(1)-SIG(2))+4.0D0*SIG(3)*SIG(3)) 
PSIG1=-0.50D0*(SIGC1)+SIG(2)  +  TEMP) 
PSIG2=-0.50D0*CSIG(1)+SIG(2)  -  TEMP) 

IF(ICODE.EQ.0)  GOTO  5002 

WRITE(6 ,978)  -EPSCl , NS) , -SIGCl) , -EPS(2, NS) , -SIGC2) , 

X  -EPS(3,NS),-SIGC3) 

WRITE(36,977) 

WRITEC36,978)  -EPS(1 , NS) , -SIGCD , -EPS(2 , NS) , -SIGC2) , 

X  -EPS(3,NS),-SIG(3) 

WRITE(24,990)  NS,  NPART,  XCELL,  YCELL 
GOTO  825 

5002  CONTINUE 

WRITE(6 ,978)  -EP(1) , -SIG(l) , -EP(2) , -SIG(2) , -EPC3) , -SIGC3) 
WRITE(36,977) 

WRITEC36,978)  -EP(1) , -SIG(l) , -EP(2) ,-SIGC2) , -EP(3) , -SIG(3) 
WRITE(24,990)  NS,  NPART,  XCELL,  YCELL 
DO  850  1=1,  NPART 

WRITEC24,994)  X(1,I),X(2,I),XC3,I) 

850  CONHNUE 

WRITE(24,991)  0.0D0,  0.0D0,  X(3,I) 

WRITE(24,991)  XA(1),  XA(2),  XC3,I) 

WRITE(24,991)  XA(1)+XB(1),  XA(2)+XBC2),  XC3,I) 

WRITE(24,991)  XB(1),  XB(2),  X(3,I) 

WRITEC24,991)  0.0D0,  0.0D0,  X(3,I) 

825  CONTINUE 

IF(NS.EQ.l)  WRITE(48,977) 

IF(ICODE.EQ.0)  WRITEC48,976)  -EP(l),-SIGCl),-EP(2),-SIGC2), 

X  -EPC3),-SIG(3) 

X  ,PSIG1,PSIG2 

IFCICODE.EQ.2)  WRITEC48 ,976)  -EPS(1,NS),-SIG(1),-EPS(2,NS), 

X  -SIG(2),-EPS(3,NS),-SIG(3) 

X  ,PSIG1,PSIG2 

1000  CONTINUE 

IF(ICODE.EQ.0)  RETURN 


STRESS  CONTROL 


c  ============ 

c 

c  sigl=sig(K-l)  sig2=sig(K) 

c  eps(i,l)=epsCK-l)  epsCi,2)=eps(K) 

7000  CONHNUE 

WRITE(6,983) 

WRITEC36,983) 

C 

C....  INPUT  MACROSCOPIC  STRESS  INCREMENTS  (DELSIG) 
C 

NSTEP  =  900 

CALL  CLEAR(DELSIG,3*NSTEP) 

CALL  CLEAR(EP,3) 

DO  5  1=1,  NSTEP 

if(i.eq.l  .or.  i.eq.2)  then 
DELSIG(1,I)  =  -0.9810d+01 
DELSIG(2,I)  =  -0.9810d+01 
DELSIG(3,I)  =  -0.00d+00 
else 

DELSIG(1,I)  =  -0.50d+00 
DELSIGC2,I)  =  0.00d+00 
DELSIGC3,I)  =  0.00d+00 
endif 
5  CONTINUE 
C 

C....  TIME  STEP  LOOP 
C 

DO  500  NS=1, NSTEP 
WRITE(6,  982)  NS 
WRITE(12,982)  NS 
WRITE(24,982)  NS 
WRITE(36,982)  NS 
C 

C....  CONVERGENCE  CRITERIA 
C 

ETOL  =  ZERO 
DO  7  1=1,  3 

ETOL  =  ETOL  +  DELSIGCI,NS)*DELSIG(I,NS) 

7  CONHNUE 

ETOL  =  1.0d-7*DSQRT(ETOL) 

C 

C....  UPDATE  CONTROL  VOLUME  CONFIGURATION 
C 

DO  25  1=1,  2 
XAOLD(I)=XA(I) 

XBOLD(I)=XB(I) 

25  CONTINUE 
C 

C....  COMPUTE  SIGbar 
C 

DO  20  1=1,  3 

FSIGCD  =  FSIGCI)  +  delsigci.ns) 


20  CONHNUE 
C 

C....  ITERAHON  LOOP 
C 

IT  =  0 

CALL  CLEAR(EPSI,3) 

450  CONTINUE 

DO  90  1=1,  2 
XACI)=XAOLD(I) 

XBCI)=XBOLDCI) 

90  CONHNUE 

IF(IT.EQ.0)  CALL  MOVECSIG,SIG2,3) 

IF(IT.EQ.0)  GOTO  499 
C 

C....  CALL  STRAIN  DRIVER  PROGRAM 
C 

WRITEC6,980) 

WRITE(36,980) 

yVRITE(6,992)  EPSI(1),EPSI(2),EPSI(3) 

WRITE(36,992)  EPSI(1),EPSI(2),EPSI(3) 

CALL  MI CROC EPSI , XA , XB , XAOLD , XBOLD , NS , NSTEP , DET, 

X  AUPPER,ALOWER,B,IDIAG,fnmax) 

WRITE(6,984) 

WRITEC36,984) 

C 

C....  UPDATE  TANGENTIAL  CONTACT  FORCES/CALCULATE  MACRO  STRESS 
C 

CALL  CLEAR(SIG,3) 

DO  35  N=1,NUMEL 
NA  =  IEN(1,N) 

NB  =  IEN(2,N) 

CALL  CNTACTCN , NA , NB , P , EK, XA , XB, DET, ITOUCH , 2 , XAOLD , 

X  XBOLD, fnmax,npl as) 

35  CONTINUE 

C 

XAxXB  =  DABS(XA(1)*XB(2)  -  XA(2)*XB(1)) 

AREA  =  XAxXB 

SIG(l)  =  SIG(1)/AREA 

SIGC2)  =  SIG(2)/AREA 

SIG(3)  =  SIG(3)/AREA 

C 

CALL  MOVE(SIGl,SIG2,3) 

CALL  MOVE(SIG2,SIG  ,3) 

C 

499  CONTINUE 

WRITEC6 ,996)  SIG(l) , SIG(2) , SIG(3) 

WRITEC36,996)  SIG(l) , SIG(2) , SIG(3) 

DO  40  1=1,  3 

DSIGCI)  =  FSIG(I)  -  SIGCI) 

40  CONHNUE 

C 

C....  COMPUTE  RESIDUAL 


c 

RESID  =  ZERO 
DO  45  1=1,  3 

RESID  =  RESID  +  DSIG(I)*DSIG(I) 

45  CONTINUE 

RESID  =  DSQRTCRESID) 

IF(IT.EQ.0)  WRITE(6,981)  ETOL 
IF(IT.EQ.0)  WRITE(12,981)  ETOL 
IF(IT.EQ.0)  WRITE(36,981)  ETOL 
WRITE(6,995)  IT+1,  RESID 
WRITE(12,995)  IT+1,  RESID 
WRITE(36,995)  IT+1,  RESID 
IF(RESID.LE.ETOL)  GOTO  700 
C 

C....  TANGENT  OPERATOR 
C 

C  TANGENT  OPTIONS  FOR  THE  START  OF  EACH  TIMESTEP: 

C  1.)  USE  THE  CONVERGED  TANGENT  FROM  THE 

C  PREVIOUS  TIMESTEP  FOR  THE  FIRST  ITERATION, 

c  IF(NS.GE.2  .AND.  IT.EQ.0)  GOTO  76 

C 

C  2.)  USE  THE  CONVERGED  TANGENT  FROM  THE 

C  PREVIOUS  TIMESTEP  FOR  THE  FIRST  TOO  ITERATIONS. 

IF(NS.GE.2  .AND.  IT.LE.l)  GOTO  76 
C 

C  3.)  UPDATE  THE  TANGENT  EVERY  OTHER  ITERATION, 

c  MDIV  =  (IT+l)/2 

c  IF(MDIV*2.EQ.(IT+1))  GOTO  76 

C 

CALL  CLEAR(SK,9) 

CALL  MOVE(SIGT,SIG2,3) 

C 

C....  DISCRETE  TANGENT  OPTIONS 

C  ITAN=0  =>  TANGENT  APPROXIMATED  BY  THE  DIFFERENCE  EXPRESSION 

C  IN  ORTEGA  AND  RHEINBOLT  (7.1.15)  PG.  185. 

C  ITAN=1  =>  TANGENT  APPROXIMATED  BY  THE  DIFFERENCE  EXPRESSION 

C  IN  ORTEGA  AND  RHEINBOLT  (7.1.16)  PG.  186. 

C  ITAN=2  =>  EXACT  TANGENT  OPERATOR 

C 

ITAN  =  0 
C 

WRITE(6,985) 

WRITE(36,985) 

DO  50  J=l,  3 
C 

IF(ITAN.EQ.l)  GOTO  8000 
C 

IF(J.GE.2)  CALL  MOVE(SIGT,SIG,3) 

IF(J.EQ.3)  CALL  MOVE(SIG,SIGl,3) 

DO  61  1=1,  3 

EPSI(I)  =  EPS(I,2) 

61  CONTINUE 


DO  63  1=1,  3 

EPSICI)  =  EPSI(I)  +  UNITY(I,K)*CEPS(K,1)-EPSCK,2)) 
63  CONHNUE 

62  CONTINUE 
GOTO  9000 
C 

8000  CONTINUE 

DO  60  1=1,  3 

EPSI(I)  =  EPS(I,2)  +  UNITYCI,J)*(EPS(J,1)-EPS(J,2)) 

60  CONTINUE 

C 

9000  CONTINUE 

write(6,*)'J  =  ',j 

writeC36,*)' J  =  ',j 

WRITEC6,992)  EPSI(l) , EPSI(2) , EPSI(3) 

WRITE(36,992)  EPSICI) , EPSI(2) , EPSI(3) 

IF(J.EQ.3  .AND.  ITAN.EQ.0)  GOTO  73 
DO  65  1=1,  2 

XA(I)  =  XAOLD(I) 

XB(I)  =  XBOLD(I) 

65  CONTINUE 

CALL  MICROCEPSI,XA,XB,XAOLD,XBOLD,NS,NSTEP,DET, 

X  AUPPER,ALOWER,B,IDIAG,fnmax) 

CALL  CLEARCSIG,3) 

DO  71  N=l,  NUMEL 
NA  =  IENC1,N) 

NB  =  IENC2,N) 

CALL  CNTACTCN , NA , NB , P , EK , XA , XB , DET , ITOUCH , 2 , XAOLD , 

X  XBOLD,fntnax,nplas) 

71  CONTINUE 

XAxXB  =  DABS(XA(1)*XB(2)  -  XAC2)*XB(1)) 

AREA  =  XAxXB 

SIGCl)  =  SIG(1)/AREA 

SIGC2)  =  SIGC2)/AREA 

SIG(3)  =  SIGC3)/AREA 

73  CONTINUE 

WRITE(6  ,996)  SIG(1),SIG(2),SIG(3) 

WRITEC36,996)  SIG(l) , SIG(2) ,SIG(3) 

TEMP  =  ONE/(EPS(J,l)-EPS(J,2)) 

DO  74  1=1,  3 

SK(I,J)  =  (SIGCD  -  SIGT(I))*TEMP 

74  CONTINUE 
50  CONTINUE 

76  continue 
WRITE(6  ,*)'SK:  ' 

WRITE(36,*)’SK:  ' 
do  77  1=1,  3 

WRITEC6  ,988)CSKCI,J),  J=l,3) 

WRITE(36,988)(SKCI,J),  J=l,3) 

77  continue 


onio  no  nnnnnoor^  nnnnr^oo 


c  writeC36,997)  dsig(l),dsig(2),dsig(3) 

WRITE(6,986) 

WRITEC36,986) 

K  =  0 

DO  105  J=l,  3 
DO  110  1=1,  J 
K  =  K  +  1 

AUPPER(K)  =  SKCI,J) 

ALOWERCK)  =  SK(J,I) 

110  CONTINUE 
105  CONTINUE 

....  FACTORIZE  TANGENT  OPERATOR 

WRITE(36,*)’ DIAGONAL  AND  RHS  BEFORE  FACTORIZATION:' 

DO  235  1=1,  3 

VI(RITE(36,987)  I,  AUPPERCIDIAG(I)),  I,  DSIGCI) 

235  CONTINUE 

CALL  NSOLVE(AUPPER, ALOWER , DSIG , IDIAG , 3 ,. TRUE .,. FALSE . ) 
DO  238  1=1,  3 

IF(DABS(AUPPERCIDIAG(I))) . LT. 1 . 0D-9) 

X  WRITEC6,979)  I,  AUPPER(IDIAG(I)) 

IF(DABS(AUPPER(IDIAG(I))) . LT. 1 . 0D-9) 

X  WRITE(12,979)  I,  AUPPER(IDIAG(I)) 

IF(DABS(AUPPER(IDIAG(I))) . LT. 1 . 0D-9) 

X  WRITE(36,979)  I,  AUPPER(IDIAG(I)) 

238  CONTINUE 

....  FORWARD  REDUCE  AND  BACK  SUBSHTUTE 

CALL  NSOLVE(AUPPER,ALOWER, DSIG, IDIAG, 3, .FALSE. , .TRUE.) 
WRITE(36,*) 

WRITE(36,*)' DIAGONAL  AND  RHS  AFTER  FACTORIZATION:' 

DO  236  1=1,  3 

WRITE(36,987)  I,  AUPPER(IDIAG(I)),  I,  DSIG(I) 

236  CONTINUE 

DO  70  1=1,  3 

EPS(I,1)  =  EPSa.2) 

EPS(I,2)  =  EPS(I,2)  +  DSIG(I) 

EPSICI)  =  EPS(I,2) 

70  CONTINUE 

WRITE(6,992)  EPSI(l) , EPSI(2) , EPSI(3) 

WRITE(36,992)  EPSI(1),EPSIC2),EPSI(3) 

WRITE(6 ,993)  DSIG(l) , DSI6(2) , DSI6(3) 

WRITE(36,993)  DSIG(l) , DSIG(2) , DSIG(3) 

IT  =  IT  +  1 
IF(IT.LE.20)  GOTO  450 

....  NO  CONVERGENCE 

WRITE(*,100) 

WRITE(36,100) 


100  FORMATC  NO  CONVERGENCE  AFTER  20  ITERATIONS') 

STOP 

700  CONTINUE 
C 

C....  CONVERGENCE 
C 

CALL  CLEAR(EPS,6) 

DO  75  1=1,3 

EP(I)  =  EPCI)  +  EPSICD 
75  CONTINUE 
C 

C....  UPDATE  TANGENHAL  CONTACT  FORCES/CALCULATE  MACRO  STRESS 
C 

fnmax=0'.0d0 
icntact  =  0 
iplstic  =  0 
CALL  CLEAR(SIG,3) 

DO  200  N=1,NUMEL 
NA  =  IEN(1,N) 

NB  =  IENC2,N) 
nplas  =  0 

CALL  CNTA CT(N , NA , NB , P , EK , XA , XB , DET, ITOUCH , 1 , XAOLD 
X  ,XBOLD,fnmax, nplas) 
icntact  =  icntact  +  itouch 
iplstic  =  iplstic  +  nplas 
200  CONTINUE 
C 

XAxXB  =  DABSCXA(1)*XB(2)  -  XA(2)*XB(1)) 

AREA  =  XAxXB 
SIGCl)  =  SIG(1)/AREA 
SI6(2)  =  SIG(2)/AREA 
SIGC3)  =  SIG(3)/AREA 
c 

write(6,*)'fnmax  =  '.fnmax 
write(24,*)'fnmax  =  ' ,fnmax 
writeC24,*)' number  of  contacts  =  '.icntact 
write(24,*)' number  of  plastic  contacts  =  '.iplstic 
C 

C....  UPDATE  COORDINATES  OF  PARTICLE  CENTROIDS 
C 

DO  80  I=1,NPART 

X(1,I)  =  X(1,I)  +  DINC(1,I) 

X(2,I)  =  X(2,I)  +  DINC(2,I) 

C 

C....  CHECK  FOR  IMAGE  PARTICLE 
C 

AA  =  (XB(2)*X(1,I)-XB(1)*X(2,I))/DET 
BB  =  (XAC1)*X(2,I)-XAC2)*XC1,I))/DET 

IF(AA. GE. ZERO. AND. AA.LE. ONE. AND. BB.GE. ZERO. AND. BB.LE. ONE) 
X  GOTO  80 
C 

C....  CONVERT  IMAGE  PARHCLE  TO  REAL 


o  r)  r>  o  o  o 


c 

IF(AA.LT.ZERO)  AA=AA+ONE 
IF(AA.GT.ONE)  AA=AA-ONE 
IF(BB.LT.ZERO)  BB=BB+ONE 
IF(BB.GT.ONE)  BB=BB-ONE 
X(1,I)  =  AA*XAC1)  +  BB*XB(1) 

X(2,I)  =  AA*XA(2)  +  BB*XB(2) 

80  CONTINUE 

....  CHECK  FOR  LOCALIZATION 
CALL  LOCALI(SK,NS) 

....  OUTPUT 

PRINCIPAL  STRESSES 

TEMP=DSQRT((SIG(1)-SIGC2))*(SIG(1)-SIG(2))+4.0D0*SIG(3)*SIG(3)) 
PSIG1=-0.50D0*(SIGC1)+SIG(2)  +  TEMP) 
PSIG2=-0.50D0*(SIGC1)+SIG(2)  -  TEMP) 

WRITEC6,978)  -EP(1) , -SIG(l) , -EPC2) , -SIG(2) , -EP(3) , -SIG(3) 
WRITE(36,977) 

WRITE(36, 978)  -EP(1) , -SIG(l) , -EPC2) , -SIGC2) , -EP(3) , -SIG(3) 
IF(NS.EQ.l)  WRITE(48,977) 

WRITE(48 , 976)  -EP(1) , -SIG(l) , -EP(2) , -SIGC2) , -EP(3) , -SIG(3) 

X  ,PSIG1,PSIG2 

WRITEC24,990)  NS,  NPART,  XCELL,  YCELL 
DO  85  1=1,  NPART 

WRITE(24,994)  xa,I),X(2,I),X(3,I) 

85  CONTINUE 

WRITE(24,991)  0.0D0,  0.0D0,  XC3,I) 

WRITEC24,991)  XA(1),  XAC2),  X(3,I) 

WRITE(24,991)  XA(1)+XB(1),  XA(2)+XB(2),  X(3,I) 

WRITEC24,991)  XB(1),  XB(2),  X(3,I) 

WRITE(24,991)  0.0D0,  0.0D0,  X(3,I) 

C 

500  CONHNUE 
C 

976  FORMAT(8E16.8) 

977  FORMATC 'MACROSCOPIC  STRESS’ ,/,5X, 

X  'EPSICD'  ,7X,  'SIGCD*  ,6X,  ’EPSI(2)'  ,7X,  ’SIGCZ)' , 

X  6X,'EPSI(3)',7X,'SIG(3)') 

978  F0RMAT(6E16.8) 

979  FORMATC'ZERO  ON  DIAGONAL—  EQ.  #',I3,'  DIAG  =  •,E20.10) 

980  FORMATC//, ’==========','  STRAIN  CONTROL', 

X  2X , ' _ ' ) 

981  FORMATC 'CONVERGENCE  CRITERION  =  ',E16.8) 

982  FORMATC//, 'TIMESTEP  ',13) 

983  FORMATC///, '===============','  STRESS  CONTROL', 

Y  9Y  • _ 

984  FORMATC  '=================================. 

985  FORMATC// ,  ’  ===== - ===== ' ,  ’  TANGENT  OPERATOR ' , 


X  2X , ' _ ' ) 

986  FORMATC  '=:1=============================== 

987  FORMATC 'DIAGC,  13,')  =  '  ,E20.10,7X, 'BC  ,13, ')  =  ',£20. 10) 

988  FORMATC3F20.10) 

990  FORMATC/, 'CONTROL  VOLUME  CONFIGURATION  AT  THE  END  OF  TIMESTEP 
,13, 

X  /,I3,2F20.14) 

991  FORMATC3F20.10) 

992  FORMATC 'STRAIN  INCREMENT: ' ,/,6X, 'EPSICl)' ,14X, 'EPSIC2)' , 

X  14X,'EPSIC3)',/,E16.8,5X,E16.8,5X,E16.8) 

993  FORMATC 'DELTACSTRAIN  INCREMENT): ' ,/,6X, 'DEPSCD' ,14X, 'DEPSC2)' 

X  14X, 'DEPSC3)’ ,/,E16.8,5X,E16.8,5X,E16.8) 

994  FORMATC3F20.10) 

995  FORMATCI2,'  STRESS  RESIDUAL  =  ',E16.8) 

996  FORMATC 'STRESS: ' ,/,6X, 'SIGCl)' ,14X, 'SIGC2)' , 

X  14X, 'SI6C3)' ,/,E16.8,5X,E16.8,5X,E16.8) 

997  FORMATC 'STRESS  RESIDUAL  VECTOR: ' ,/,6X, 'DSIGCD' ,14X, 'DSIGC2)' , 

X  14X, 'DSIGC3)' ,/,E16.8,5X,E16.8,5X,E16.8) 


or^r^nnn 


SUBROUTINE  MI CRO( EPSI , XA , XB , XAOLD , XBOLD , NS , NSTEP , DET, 

X  AUPPER,ALOWER,B,IDIAG,fnmax) 

....  STRAIN  DRIVEN  PROBLEM  DRIVER  PROGRAM 

IMPLICIT  REAL*8(A-H,0-Z) 

....  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 

DIMENSION  EPSI(3) , XA0LDC2) , XB0LD(2) 

DIMENSION  EKC6 , 6) , P(6) , XA(2) , XBC2) 

DIMENSION  AUPPER(l) ,AL0WERC1) ,B(1), lOIAG(l) ,BDPLC1000) 
COMMON  /MACRO  /  FSIG(3),SIGC3),EPC3) 

COMMON  /PARTCL/  X(3,1000),NPART,RMIN 

COMMON  /INFO  /  IENC2, 20000), 10(3, 1000), LM(6, 20000), 

X  NUMEL,NEQ,NSTIFF 

COMMON  /DSDATA/  0(3, 1000), 00(3, 1000), DINC(3, 1000), 

X  DDINC(3,1000) 

COMMON  /FDATA  /  FTN(20000),ALPHA(20000),ICON(20000) 
COMMON  /ELDATA/  SKN,SKT,A,TANPHI,HP 

DATA  ZERO , ONE , ETOL/0 . 0D0 , 1 . 0D0 , 1 . 0D-8/ 

....  INITIALIZE  Dbap 

DO  5  I=1,NPART 

DINC(1,I)  =  EPSI(1)*X(1,I)  +  EPSI(3)*X(2,I) 

DINC(2,I)  =  EPSI(3)*X(1,I)  +  EPSI(2)*X(2,I) 

DINC(3,I)  =  0.0D0 
5  CONTINUE 

TEMP  =  (ONE+EPSI(l))*XA(l)  +  EPSI(3)*XA(2) 

XA(2)  =  EPSI(3)*XA(1)  +  (ONE+EPSI(2))*XA(2) 

XA(1)  =  TEMP 

TEMP  =  (ONE+EPSI(l))*XB(l)  +  EPSI(3)*XB(2) 

XB(2)  =  EPSI(3)*XB(1)  +  (ONE+EPSI(2))*XB(2) 

XB(1)  =  TEMP 

DET  =  XA(1)*XB(2)-XB(1)*XA(2) 

....  PSEUDO  TIME  STEP  LOOP 

DT  =  1.0D0 
DO  400  NT=1,  1 
WRITE(6,980)  NT 
WRITE(12,980)  NT 
WRITE(24,980)  NT 
CALL  CLEAR(DDINC,3*NPART) 

....  NEWTON  ITERATION  LOOP 


ITER  =  0 


10  CONnNUE 

ITER  =  ITER  +  1 
c  1NRITE(6,996)  ITER 

c  WRITE(12,996)  ITER 

C 

701  CONTINUE 
ALPHK  =  ONE 

C 

702  CONHNUE 

CALL  CLEAR(AUPPER,NSTIFF) 

CALL  CLEARCALOWER.NSTIFF) 

CALL  CLEAR(B,NEQ) 

C 

C....  FORM  RHS  VECTOR  AND  TANGENT  OPERATOR 
C 

iflag  =  0 
DO  15  N=1,NUMEL 
CALL  CLEAR(P,6) 

CALL  CLEAR(EK,36) 

NA  =  IEN(1,N) 

NB  =  IEN(2,N) 

CALL  CNTACT(N , NA , NB , P , EK, XA , XB, DET, ITOUCH , ITASK , XAOLD , 

X  XBOLD,fntnax,nplas) 
c  If(itouch.eq.l)  then 
c  if(iflag.eq.l)  goto  104 
c  iflag  =  1 

c  write(12,*)  'n  =  ',n 

c  do  303  ii=l,6 

c  write(12,223)  (ek(ii,jj),jj=l,6),  p(ii) 
c  223  formatC7fl0.3) 
c  222  foriiiatC6fl0.3) 
c  303  continue 
c  endif 
c  104  continue 

IF(ITOUCH.EQ.0)  GOTO  15 

CALL  ADDSTF(AUPPER,B,ALOWER,EK,P,IDIAG,LMCl,N)) 

15  CONTINUE 
C 

C....  PARABOLIC  REGULARIZATION 
C 

IF(ITER.EQ.l)  WRITE(6, *)’***  NO  PSEUDO  TIMESTEPPING  ***' 
GOTO  35 

IFCITER.EQ.I)  goto  35 
DO  30  N=l,  NPART 

CALL  PAREG(N , AUPPER , B , ALOWER, IDIAG , DT) 

30  CONTINUE 
35  CONTINUE 
C 

C....  COMPUTE  NORM  OF  RESIDUAL  AND  CHECK  CONVERGENCE 
C 

RESIDl  =  DSQRT(DOT(B,B,NEQ)) 

IF(ITER.EQ.l)  TOLER=RESIDl*ETOL 


nnnn 


_  STEPLENGTH  DETERMINATION/LINE  SEARCH 

FOR  NO  LINE  SEARCH  GOTO  705 
go  to  705 
706  continue 

IFCRESIDl . LT. RESID0 .OR. ITER . EQ. 1)  THEN 
IFCITER.EQ.l)  GOTO  705 
WRITE(12,989) 

DO  22  J=l,  NPART 

WRITE(12,985)  J,  (DDINC(I,J),  1=1,3),  (DINC(I, J),I=1,3) 

22  .  CONHNUE 
GOTO  705 
ENDIF 

ALPHK  =  0.50D0*ALPHK*ALPHK*RESID5/CRESID1+RESID5*ALPHK-RESID0) 
WRITE(*,777)  ALPHK 
WRITE(36,777)  ALPHK 
777  FORMATC  ALPHK  =  •,F20.10) 

DO  17  1=1,3 
DO  16  J=l, NPART 
K  =  IDCI,J) 

IF(K.GT.0)  DINC(I,J)=DCI,J)-ALPHK*BDPL(K) 

IF(K. GT. 0)  DDINCCI , J)=DD(I , J)-ALPHK*BDPL(K) 

16  CONHNUE 

17  CONHNUE 
GO  TO  702 

705  CONTINUE 

IFCITER.EQ.l)  WRITE(6,981)  TOLER 
IFCITER.EQ.l)  WRITEC36,981)  TOLER 
WRITEC6,995)  ITER,  RESIDl 
Vi/RITEC36,995)  ITER,  RESIDl 

IFCRESIDl. LT.TOLER.OR. RESIDl. LT.1.0D-7)  GOTO  400 
RESID0  =  RESIDl 
RESID5  =  RESID0 

....  FACTORIZE  TANGENT  OPERATOR 

c  WRITEC12,*)' DIAGONAL  AND  RHS  BEFORE  FACTORIZATION:' 
c  DO  235  1=1,  NEQ 

c  WRITEC12,987)  I,  AUPPERCIDIAGCD) ,  I,  BCI) 

c  235  CONTINUE 

CALL  NSOLVECAUPPER,ALOWER,B,IDIA6, NEQ,. TRUE. , .FALSE.) 

DO  238  1=1,  NEQ 

c  IFCDABSCAUPPERCIDIAGCI))) . LT. 1. 0D-9) 

c  X  WRITEC6,979)  I,  AUPPERCIDIAGCI)) 

c  IFCDABSCAUPPERCIDIAGCI))) . LT. 1 . 0D-9) 

c  X  WRITEC12,979)  I,  AUPPERCIDIAGCI)) 

IFCDABSCAUPPERCIDIAGCI))) . LT. 1 . 0D-9) 

X  WRITEC36,979)  I,  AUPPERCIDIAGCI)) 

238  CONTINUE 
C 

C....  FORWARD  REDUCE  AND  BACK  SUBSTITUTE 


CALL  NSOLVE(AUPPER,ALOWER,B,IDIAG,NEQ,. FALSE.,. TRUE.) 
c  WRITEC12,*) 

c  WRITE(12,*)' DIAGONAL  AND  RHS  AFTER  FACTORIZATION:' 

c  DO  236  1=1,  NEQ 

c  WRITE(12,987)  I,  AUPPERCIDIAG(I)),  I,  B(I) 

236  CONTINUE 

....  INTERMEDIATE  UPDATE  FOR  DISPLACEMENT  INCREMENT 

CALL  M0VECD,DINC,3*NPART) 

CALL  M0VE(DD,DDINC,3*NPART) 

CALL  MOVE(BDPL,B,NEQ) 

DO  25  1=1,3 
DO  20  J=1,NPART 
K  =  IDCI,J) 

IF(K.GT.0)  DINC(I,J)=DINC(I,J)-B(K) 

IF(K.GT.0)  DDINC(I,J)=DDINC(I,J)-BCK) 

20  CONHNUE 
25  CONTINUE 

IF  (ITER.LT.45)  GOTO  10 

....  NO  CONVERGENCE 

WRITEC*,100) 

WRITEC36,100) 

100  FORMATC  NO  CONVERGENCE  AFTER  45  ITERATIONS.’) 

STOP 

400  CONTINUE 
C 

979  FORMATC'ZERO  ON  DIAGONAL—  EQ.  #’,I3,'  DIAG  =  ',E20.10) 

980  FORMATC 'PSEUDO  TIMESTEP  ’,13) 

981  FORMATC 'CONVERGENCE  CRITERION  =  ',E16.8) 

985  FORMATCI3,2X,3E12.4,’  l',3E12.4) 

987  FORMATC 'DIAGC,  13,')  =  ’  ,E20.10,7X, 'BC  ,13, ')  =  ’,E20.10) 

989  FORMATC/, 'PRT.#' ,2X, ’DCDINCCl))’ ,2X, 'DCDINCC2))' ,2X, ’DCDINCC3))' , 
X  IX , ' I ' ,4X , ' DINCCl) ' , 5X , ’ DINCC2) ' , 5X , ' DINCC  3) ' ) 

995  FORMATCI2,’  RESIDUAL  =  ',E16.8) 

996  FORMATC//, 'ITERATION  =  ’,12) 

C 

RETURN 

END 


S UBROUTI N E  CNTA  Cr(N , NA , NB , P , EK , XA , XB , DET , ITOUCH , ITASK , 

X  XAOLD,XBOLD,fnmax,nplas) 

PROGRAM  TO  FORM  ELEMENT  FORCE  VECTOR  AND  SHFFNESS  MATRIX 

IMPLICIT  REAL*8(A-H,0-Z) 

. . . .  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 

DIMENSION  P(l) , EK(6 , 1) ,XA(1) , XB(1) , XAOLDCl) , XBOLD(l) 
DIMENSION  XN1(2 , 2) , XL0C2) , XL1(2) , XNB(2) 

DIMENSION  C(2,2),B(2,4),CB(2,4) 

COMMON  /MACRO  /  FSIG(3),SIG(3),EPC3) 

COMMON  /PARTCL/  X(3,1000),NPART,RMIN 

COMMON  /DSDATA/  D(3, 1000), DD(3, 1000). DINCC3, 1000), 

X  DDINC(3,1000) 

COMMON  /FDATA  /  FTN(20000),ALPHA(20000),ICON(20000) 

COMMON  /ELDATA/  SKN,SKT,A,TANPHI,HP 
DATA  ZERO , HALF , ONE/0 . 0D0 , 0 . 50D0 , 1 . 0D0/ 

....  ITASK  OPnONS: 

ITASK=0  =>  CALCULATE  ELEMENT  CONTRIBUTION  TO  MICRO  STIFFNESS 
ITASK=1  =>  UPDATE  MACRO  STRESS  AND  CONTACT  FORCES 
ITASK=2  =>  CALCULATE  MACRO  STRESS  WITHOUT  ANY  UPDATES 

....  lOPT  OPnONS: 

IOPT=0  =>  ALL  PARTICLE  IN  CONTACT  WITH  EACH  OTHER  AT  Tn 
IOPT=l  =>  NEW  CONTACTS  FORM  AT  Tn+1. 

lOPT  =  0 

RA  =  X(3,NA)/100.0D0 
RB  =  X(3,NB)/100.0D0 
RPLUSR  =  (RA  +  RB)/100.0D0 
RA  =  X(3,NA) 

RB  =  X(3,NB) 

RPLUSR  =  (RA  +  RB) 

XN1(1,1)  =  X(1,NA)  +  DINC(1,NA) 

XN1(2,1)  =  X(2,NA)  +  DINC(2,NA) 

XN1(1,2)  =  X(1,NB)  +  DINC(1,NB) 

XN1(2,2)  =  X(2,NB)  +  DINC(2,NB) 

THETAA  =  DINC(3,NA) 

THETAB  =  DINC(3,NB) 

. . . .  CHECK  REAL  AND  IMAGE  PARTICLES 

XLl(l)  =  XN1(1,2)  -  XN1(1,1) 

XL1(2)  =  XN1(2,2)  -  XN1(2,1) 

AA  =  (XB(1)*XL1(2)-XB(2)*XL1(1))/DET 


ooo  ooo  oor^  ooo 


•  BB  =  CXAC2)*XL1C1)-XA(1)*XL1(2))/DET 

SIGNA  =  ZERO 

IF(AA.NE.ZERO)  SIGNA=AA/DABSCAA) 

SIGNB  =  ZERO 

IFCBB.NE.ZERO)  SIGNB=BB/DABSCBB) 

C 

^  ITOUCH  =  0 

•  DO  40  1=0,1 

DO  35  J=0,1 
DI  =  DREAL(I) 

DJ  =  DREAL(J) 

XNB(l)  =  XN1(1,2)  +  DI*SIGNA*XAC1)  +  DJ*SIGNB*XB(1) 
XNB(2)  =  XN1(2,2)  +  DI*SIGNA*XA(2)  +  DJ*SIGNB*XBC2) 

m  c 

XL0(1)  =  XCl.NB)  +  DI*SIGNA*XA0LDC1)  +  DJ*SIGNB*XB0LD(1) 
XL0(2)  =  X(2,NB)  +  DI*SIGNA*XA0LD(2)  +  DJ*SIGNB*XB0LD(2) 

....  BRANCH  VECTORS 

XLl(l)  =  XNB(l)  -  XN1(1,1) 

XL1(2)  =  XNB(2)  -  XN1C2,1) 

XL0(1)  =  XL0C1)  -  X(1,NA) 

XL0(2)  =  XL0(2)  -  X(2,NA) 

DELTA0  =  DSQRTCDOT(XL0,XL0,2)) 

DELTAl  =  DSQRTCD0T(XL1,XL1,2)) 

....  CHECK  CONTACT 

DELTA  =  RPLUSR  -  DELTAl 
IF(DELTA.GT.ZERO)  GOTO  50 
35  CONTINUE 
40  CONTINUE 

....  NO  CONTACT 

IF(ITASK.GE.0)  RETURN 
ICON(N)  =  0 
FTN(N)  =  ZERO 
ALPHACN)  =  A 
RETURN 


....  CONTACT 


50 


CONTINUE 

IF(NA.EQ.7  .OR.  NB.EQ.7)  WRITEC60,*)'  NA  =  ',NA,'  NB  =  ',NB, DELTA 
WPite(36,*)' CONTACT  BETWEEN  PARTICLES  ',NA,'  AND  ’ ,NB,  DELTA 
ITOUCH  =  1 
EPSl  =  0.1D0*RMIN 
HEPS  =  ONE 
DHEPS  =  ZERO 

IFCDELTA. LT. EPSl)  HEPS=DELTA/EPS1 
IF(DELTA.LT.EPSl)  DHEPS=ONE/EPSl 


. . . .  CONTACT  STIFFNESS 

PI=DATAN(1.0D0)*4.0[)0 
E  CN/ctnA2) 

E  =  5.00d4 

NU  =  0.20D0 

R  =  1. 000/(1. 0D0/RA  +  1.0D0/RB) 

Po  =  3.0d0 

TEMP  =  RA*RB*PI*E/(4.0D0*R*(1.0D0-NU*NU)*Po) 

SKN  =  PI*E/(4.0D0*C1.0D0-NU*NU)*CDLOG(TEMP)+1.0D0)) 

SKT  =  SKN 


. . . .  NORMALIZE  BRANCH  VECTORS  # 

DO  55  1=1,2 

XL0(I)  =  XL0(I)/DELTA0 
XLICD  =  XL1CI)/DELTA1 
55  CONTINUE 

. . . .  ROTATION  ANGLE  • 

THETAC  =  DASINCXL0C1)*XL1(2)  -  XL0C2)*XL1C1)) 

WRITE(12,*)' THETA  C  =  ',  THETAC 

. . . .  NORMAL  CONTACT  FORCE 

FN  =  -HEPS*SKN*DELTA 

. . . .  TANGENTIAL  CONTACT  FORCE 

GAMMA  =  RA*CTHETAC-THETAA)  +  RB*CTHETAC-THETAB) 

FTRIAL  =  FTN(N)  +  HEPS*SKT*GAMMA 

YIELD  =  ALPHA(N)  -  FN*TANPHI  ♦ 

WRITEC 12,*) 'FTRIAL  =  ', FTRIAL 
WRITE(12, *)' YIELD  =  ’, YIELD 
IFCIOPT.EQ.l  .AND.ICON(N).EQ.0)  THEN 
GPLAS  =  ZERO 
FT  =  FTRIAL 

C(l,l)  =  -SKN*CHEPS+DELTA*DHEPS)  a 

C(l,2)  =  ZERO 
C(2,l)  =  ZERO 
CC2,2)  =  ZERO 
GO  TO  90 
ENDIF 

IF(DABSCFTRIAL).GT. YIELD)  GOTO  70  ^ 

....  ELASTIC  PROCESS;  C  =  LOCAL  STRESS-STRAIN  MATRIX 


c 

c 


WRITE(12,*)' ELASTIC  PROCESS' 
writeC*,*)  'elastic  process' 
nplas  =  0 


mr^r^nonnn  nn  ooo 


GPLAS  =  ZERO 
FT  =  FTRIAL 

CC1,1)  =  -SKN*(HEPS+DELTA*DHEPS) 

CC1,2)  =  ZERO 
0(2,1)  =  SKT*GAMMA*DHEPS 
CC2,2)  =  HEPS*SKT 
c  WRITEC 12,*) 'GAMMA  =  '.GAMMA 

c  WRITE(12,999)N,NA,NB,FT,FN 

GO  TO  90 

....  PLASTIC  PROCESS;  C  =  LOCAL  STRESS-STRAIN  MATRIX 

70  CONTINUE 
nplas  ='  1 

if(itask.eq.l)  WRITEC24,*)' PLASTIC  PROCESS' 
write(*,*)  'plastic  process' 

SIGN  =  FTRIAL/DABS( FTRIAL) 

GPLAS  =  SIGN*CDABSCFTRIAL)-YIELD)/(HEPS*SKT+HP) 

FT  =  FTRIAL  -  HEPS*SKT*GPLAS 
TEMP  =  ONE  +  HP/(SKT*HEPS) 

C(l,l)  =  -SKN*(HEPS+DELTA*DHEPS) 

CC1,2)  =  ZERO 

C(2,l)  =  HP*DHEPS*(GAMMA-GPLAS)/HEPS 
X  +SIGN*SKN*TANPHI*(HEPS+DELTA*DHEPS) 

CC2,1)  =  C(2,1)/TEMP 

C(2,2)  =  HEPS*SKT*(ONE  -  ONE/TEMP) 

WRITEC12,*)'SIGN  =  ' .SIGN 
WRITEC 12,*) 'GAMMA  =  '.GAMMA 
WRITEC 12,*) 'PLASTIC  CAMMA  =  '.gplas 
WRITEC12,*)' ELASTIC  GAMMA  =  ' ,gamma-gplas 
WRITEC12 ,999)N , NA , NB , FT, FN 

....  UPDATE  VARIABLES 

90  CONHNUE 

IFCITASK.EQ.0)  GOTO  100 
c  writeC*,102) 

102  formate'  i  am  here  at  update...') 

IFCIOPT.EQ.l  .AND.  ICON(N). EQ.0)  FT  =  ZERO 
FX  =  XL1(1)*FN  -  XL1(2)*FT 
FY  =  XL1C2)*FN  +  XL1C1)*FT 
C 

SIG(l)  =  SIGCD  +  FX*XL1(1)*DELTA1 
SIG(2)  =  SIG(2)  +  FY*XL1(2)*DELTA1 
SIG(3)  =  SIGC3)  +  0.50D0*DELTA1*CFX*XL1C2)+FY*XL1C1)) 
IF(ITASK.EQ.l)  THEN 
ICONCN)  =  1 
FTNCN)  =  FT 

ALPHA(N)  =  ALPHACN)+HP*SIGN*GPLAS 
c 

c  WRITE(24,505)  NA,NB,FX,FY,FT,FN 

c  R=1.0D0/(1.0D0/RA  +  1.0D0/RB) 


c  e=2.0d0*E/(1.0d0-nu*nu) 

c  a=8 . 0*f n*  r/3 . 14160d0/e 

c  a=a*a/r 

c  write(24,*)'a/R  =  a 

if(dabs(fn) . gt.fnmax)fnmax=dabs(fn) 

RETURN 

ENDIF 

IF(ITASK.EQ.2)  RETURN 

505  FORMATC  ELEM. '  ,13, ,13, '  FX=’,F10.5, 

X  '  FY=',F10.5,'  FT=',E12.4,'  FN=',E12.4) 

RETURN 
C 

C....  ASSEMBLE  INTERNAL  FORCE  VECTOR 
C 

100  CONHNUE 

P(l)  =  XL1(1)*FN  -  XL1(2)*FT 
P(2)  =  XL1C2)*FN  +  XL1(1)*FT 
PC3)  =  RA*FT 
PC4)  =  -P(l) 

P(5)  =  -P(2) 

PC6)  =  RB*FT 
C 

C....  B  =  STRAIN/DISPLACEMENT  MATRIX 
C 

BC1,1)  =  XLl(l) 

B(l,2)  =  XL1(2) 

B(l,3)  =  ZERO 
B(l,4)  =  ZERO 

TANC  =  DTAN(THETAC)/DELTA1 
SECC  =  ONE/(DCOS(THETAC)*DELTAl) 

BC2,1)  =  (XL1(1)*TANC  +  XL0C2)*SECC)*RPLUSR 
B(2,2)  =  (XL1C2)*TANC  -  XL0C1)*SECC)*RPLUSR 
BC2,3)  =  -RA 
B(2,4)  =  -RB 
c  WRITE(12,*)'B' 
c  DO  112  1=1,  2 
c  WRITE(12,994)  (BCI,J),  J=l,  4) 

c  112  CONTINUE 
C 

C....  CB  =  DISPLACEMENT  GRADIENT  OF  LOCAL  FORCE  VECTOR 
C 

DO  110  1=1,2 
DO  110  J=l,4 

CB(I,J)  =  Ca,l)*B(l,J)  +  Ca.2)*BC2,J) 

110  CONHNUE 
c  WRITEC12,*)’CB' 
c  DO  111  1=1,  2 

c  WRITE(12,994)  (CB(I,J),  J=l,  4) 

c  111  CONTINUE 
C 

C....  DISPLACEMENT  GRADIENT  OF  UNIT  NORMAL 
C 


c 

c.. 

c 


c 

c 


120 

C 

C.... 

C 


C 


130 


C 


C 

c 

c 

c 

c  73 
c 
c 
c 

c  74 


SLll  =  (XL1(1)*XL1(1)-0NE)/DELTA1 
SL12  =  XL1(2)*XL1(1)/DELTA1 
SL22  =  (XL1(2)*XL1(2)-0NE)/DELTA1 

C  =  DISPLACEMENT  GRADIENT  OF  GLOBAL  FORCE  VECTOR 

C(l,l)  =  XL1(1)*CBC1,1)  -  XL1(2)*CBC2,1)  +  SL11*FN  -  SL12*FT 

C(l,2)  =  XL1(1)*CB(1,2)  -  XL1C2)*CB(2,2)  +  SL12*FN  -  SL22*FT 

C(2,l)  =  XL1(2)*CBC1,1)  +  XL1(1)*CB(2,1)  +  SL12*FN  +  SL11*FT 

C(2,2)  =  XL1(2)*CB(1,2)  +  XL1(1)*CB(2,2)  +  SL22*FN  +  SL12*FT 

ELEMENT  STIFFNESS  MATRIX 

DO  120  1=1,2 
DO  120  J=l,2 

EK(I  ,J  )  =  EKCI  ,J  )  +  C(I,J) 

EK(I+3,J  )  =  EK(I+3,J  )  -  C(I,J) 

EK(I  ,J+3)  =  EK(I  ,J+3)  -  C(I,J) 

EK(I+3,J+3)  =  EK(I+3,J+3)  +  C(I,J) 

CONTINUE 

C  =  ROTATION  GRADIENT  OF  GLOBAL  FORCE  VECTOR 

CC1,1)  =  XL1(1)*CB(1,3)  -  XL1(2)*CBC2,3) 

C(l,2)  =  XL1(1)*CB(1,4)  -  XL1C2)*CB(2,4) 

C(2,l)  =  XL1(2)*CB(1,3)  +  XL1(1)*CB(2,3) 

CC2,2)  =  XL1(2)*CB(1,4)  +  XL1(1)*CBC2,4) 

DO  130  1=1,2 

EK(I  ,3)  =  EK(I  ,3)  +  C(I,1) 

EK(I+3,3)  =  EKCI+3,3)  -  0(1,1) 

EK(I  ,6)  =  EK(I  ,6)  +  C(I,2) 

EKCI+3,6)  =  EKCI+3,6)  -  C(I,2) 

EKC3,I  )  =  EKC3,I  )  +  RA*CBC2,I) 

EK(3,I+3)  =  EK(3,I+3)  -  RA*CB(2,I) 

EK(6,I  )  =  EK(6,I  )  +  RB*CB(2,I) 

EK(6,I+3)  =  EK(6,I+3)  -  RB*CB(2,I) 

CONTINUE 

EK(3,3)  =  EK(3,3)  +  RA*CB(2,3) 

EK(3,6)  =  EKC3,6)  +  RA*CBC2,4) 

EK(6,3)  =  EKC6,3)  +  RB*CB(2,3) 

EK(6,6)  =  EK(6,6)  +  RB*CB(2,4) 

WRITEC12,*)' STIFFNESS  MATRIX' 

DO  73  1=1,  6 

WRITE(12,996)(EKCI,J),  1=1,  6) 

CONTINUE 

WRITEC12,*)' ELEMENT  INTERNAL  FORCE  VECTOR’ 

DO  74  1=1,  6 

WRITE(12,995)  PCI) 

CONTINUE 


c 

994  F0RMATC4F14.6) 

995  F0RMAT(F12.4) 

996  F0RMATC6F12.4) 

999  FORMATC ’ ELEMENT  NUMBER  ',15,'  CONNECTING  PARTICLES  ',13, 
X  '  AND  ',13,/, 'TANGENTIAL  FORCE  =  ',E16.8,/, 

X  'NORMAL  FORCE  =  ',E16,8,/) 


RETURN 

END 


r^oor^or>nnnnn  000000000  o  000  000 


c 


SUBROUTINE  LOCALI(C,NS) 


. . . .  INVERSE  PROBLEM  DRIVER  PROGRAM 
IMPLICIT  REAL*8(A-H,0-Z) 

. . . .  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 

DIMENSION  C(3,3),X(3) 

DATA  ZERO, ONE/0. 0D0,1.0D0/ 

HALF  =  1.0D0/2.0D0 
THIRD  =  1.0D0/3.0D0 
RAD  =  4.0D0*DATANC1.0D0)/180.0D0 

_  CALCULATE  THE  DETERMINATE  OF  THE  ACOUSTIC  TENSOR. 

A0,  Al,  A2,  A3,  AND  A4  ARE  THE  COEFFICIENTS  OF  THE 
QUARTIC  EXPRESSION  DEFINING  THE  LOCALIZATION  CONDITION. 

F(X)  =  A4  X**4  +  A3  X**3  +  A2  X**2  +  Al  X  +  A0  =  0 
WHERE  X=TANCTHETA),  THETA  DEFINES  THE  ORIENTATION  OF 
THE  SHEAR  BAND. 

A0  =  CC1,1)*C(3,3)  -  C(1,3)*C(3,1) 

Al  =  CC1,1)*C(3,2)  +  C(1,1)*C(2,3) 

X  -  CC1,3)*CC2,1)  -  CC1,2)*C(3,1) 

A2  =  CC1,1)*C(2,2)  +  C(1,3)*C(3,2)  +  C(3,1)*C(2,3) 

X  -  CC1,2)*CC3,3)  -  C(1,2)*C(2,1)  -  C(3,3)*C(2,1) 

A3  =  C(1,3)*C(2,2)  +  C(3,1)*C(2,2) 

X  -  C(1,2)*CC2,3)  -  C(3,2)*C(2,1) 

A4  =  C(3,3)*C(2,2)  -  C(2,3)*C(3,2) 

a0  =  1.0d0 
al  =  -240.0334830d0 
aZ  =  1.0d0 
a3  =  -10.0d0 
a4  =  1.0d0 

....  FIND  THE  MINIMA  OF  THE  LOCALIZATION  CONDIHON.  THE 
MINIMA  OCCURS  AT  THE  ROOTS  OF  THE  DERIVATIVE  OF  THE 
LOCALIZATION  CONDITION  —  A  CUBIC  EQUATION. 

F'(X)  =  X**3  +  P  X**2  +  Q  X  +  R  =  0 

P  =  A3/A4*3.0D0/4.0D0 
Q  =  A2/A4/2.0D0 
R  =  A1/A4/4.0D0 

....  FIND  THE  ROOTS  OF  F'(X). 

A  =  (3.0D0*Q  -  P*P)*THIRD 

B  =  (2.0D0*P*P*P  -  9.0D0*P*Q  +  27.0D0*R)/27.0D0 


QQ  =  B*B/4.0D0  +  A*A*A/27.0D0 
IF(DABSCQ®.LT.1.0D-8)  QQ  =  ZERO 
writeCe,*)  'QQ  =  ’,qq 

. . . .  ANALYSIS  OPTIONS 

QQ>0  =>  ONE  REAL  ROOT  AND  TOO  CONJUGATE  COMPLEX  ROOTS; 

QQ=0  =>  THREE  REAL  ROOTS  OF  WHICH  TOO  AT  LEAST  TOO  ARE  EQUAL; 
QQ<0  =>  THREE  REAL  AND  UNEQUAL  ROOTS. 

....  NOTE  THE  CHANGE  IN  VARIABLES  Y  =  X  -  P/3. 

TEMP  =  P*THIRD 

....  COMPUTE  ONLY  THE  REAL  PARTS  OF  ALL  THE  ROOTS.  NOTE  THAT 
THE  COMPLEX  ROOTS  WILL  NOT  BE  USED  IN  THE  LOCALIZATION 
ANALYSIS  BUT  THEIR  REAL  PARTS  WILL  BE  CALCULATED  ANYWAYS. 

. . . .  FOR  QQ>0  OR  QQ=0 

IF(QQ.LT.ZERO)  GOTO  500 

. . . .  FOR  QQ<0 
TEMP2  =  ONE 

TEMP3  =  -HALF*B  +  DSQRT(QQ) 

IF(TEMP3.LT.ZERO)  TEMP2  =  -ONE 
AA  =  TEMP2*(DABS(TEMP3))**THIRD 
TEMP2  =  ONE 

TEMP3  =  -HALF*B  -  DSQRT(QQ) 

IF(TEMP3.LT.ZERO)  TEMP2  =  -ONE 
BB  =  TEMP2*(DABS(TEMP3))**THIRD 
X(l)  =  AA  +  BB  -  TEMP 
X(2)  =  -(AA  +  BB)/2.0D0  -  TEMP 
X(3)  =  X(2) 

GOTO  1000 
500  CONTINUE 

THETA  =  DACOS(-HALF*B/DSQRT(DABS(-A*THIRD*A*THIRD*A*THIRD))) 

TEMP2  =  2.0D0*DSQRT(DABS(-THIRD*A)) 

X(l)  =  TEMP2*DCOS(THETA*THIRD)  -  TEMP 

XC2)  =-TEMP2*DCOS(THETA*THIRD  +  60.0D0*RAD)  -  TEMP 

X(3)  =-TEMP2*DCOS(THETA*THIRD  -  60.0D0*RAD)  -  TEMP 

1000  CONTINUE 

_  FIND  THE  MINIMIA  OF  THE  LOCALIZATION  FUNCTION  USING 

ONLY  THE  REAL  ROOTS  OF  THE  CUBIC  EQUATION. 

FMIN  =  1.0D50 
N  =  3 

IF(DABS(Q0.LT.1.0D-8)  N  =  2 
IF(QQ.GT.ZERO)  N  =  1 
DO  5  1=1,  N 

F  =  A4*X(I)*X(I)*X(I)*XCI)  +  A3*X(I)*X(I)*XCI) 

+  A2*X(I)*X(I)  +  A1*X(I)  +  A0 


X 


o  o 


IFCF.GT.FMIN)  GOTO  5 
FMIN  =  F 
XMIN  =  X(I) 

CONHNUE 


....  OUTPUT 

I FCFMIN . LT. 1 . 0D-8)  WRITE(36 , *) ' LOCALIZATION ’ 

IF(FMIN.LT.1.0D-8)  WRITE(  6,*)’ LOCALIZATION’ 

WRITE(36,990)N,  FMIN,  XMIN 
c  WRITEC  6,990)N,  FMIN,  XMIN 
WRITEC  6,993)  N 
WRITE(36,*)A4,  A3,  A2,  Al,  A0 

•  IF(NS.EQ.l)  WRITE(60,991) 

WRITE(60,992)  NS,  FMIN,  N 
WRITE(36,*)  'ROOTS:  ' 

WRITE(36,*)  XCD,  XC2),  X(3) 

990  FORMATC'NUM.  OF  REAL  ROOTS  =  ',11,/, 

X  'LOCALIZATION  FUNCTION  MINIMA  =  ',E16.8,/, 

_  X  ’TAN(THETA)  AT  MINIMA  =  ',E16.8) 

•  991  FORMATC ' LOCALIZATION  FUNC.  FOR  EACH  TIMESTEP  ', 

X  'AND  NUM.  OF  REAL  ROOTS: ',/, 'TIMESTEP' ,3X, ' LOCAL.  FUNCT. ',3X, 
X  'NUM.  REAL  ROOTS’) 

992  FORMATCI3,5X,E20.10,5X,I2) 

993  FORMATC'NUM.  OF  REAL  ROOTS  =  ’,11) 

C 

•  RETURN 

END 


ni  o  o  o  o 


c============«============================= 

SUBROUTINE  PAREG(N , AUPPER , B , ALOWER , IDIAG , DT) 

PROGRAM  TO  FORM  PARTICLE  FORCE  VECTOR  AND  STIFFNESS  MATRIX 
RESULTING  FROM  THE  DRAG/VISCOUS  FORCE. 

IMPLICIT  REAL*8(A-H,0-Z) 

....  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 

DIMENSION  AUPPER(l) , B(l) , ALOWER(l) , IDIAG(l) 

COMMON  /PARTCL/  X(3,1000),NPART,RMIN 

COMMON  /INFO  /  IEN(2, 20000), 10(3, 1000), I-MC6, 20000), 

X  NUMEL,NEQ,NSnFF 

COMMON  /DSDATA/  D(3, 1000), DD(3, 1000), DINC(3, 1000), 

X  DDINC(3,1000) 

DATA  ZERO,T1IVO/0.0D0,2.0D0/ 

C 

c  DAMP  =-1.0D2 

DAMP  =  0.0D0 
R  =  X(3,N) 

C 

DMAG  =  DDINC(1,N)*DDINC(1,N)  +  DDINC(2,N)*DDINCC2,N) 

DRT  =  DSQRT(DMAG) 
c  wpite(12,*)  n 
c  write(12,997)  drt 

IF(DRT.LT.1.0D-8)  GOTO  300 
C 

RDAMP  =  R*DAMP*DRT/(DT*DT) 

C 

DO  200  1=1,  2 
K  =  ID(I,N) 

IF(K.EQ.0)  GOTO  200 
Vise  =  RDAMP*DDINCCI,N) 

B(K)  =  B(K)  +  Vise 
c  write(12,998)  i,  vise 

200  CONTINUE 
C 

RDAMP  =  RDAMP/(DRT*DRT) 

K  =  ID(1,N) 

IF(K.EQ.0)  GOTO  210 
J  =  IDIAG(K) 

STF  =  RDAMP*CDDINC(1,N)*DDINC(1,N)  +  DMAG) 
c  write(12,996)  1,  1,  stf 

AUPPER(J)  =  AUPPER(J)  +  STF 
ALOWER(J)  =  ALOWERCJ)  +  STF 
C 

210  CONTINUE 
K  =  IDC2,N) 

IF(K.EQ.0)  GOTO  220 
J  =  IDIAG(K) 

STF  =  RDAMP*CDDINC(2,N)*DDINC(2,N)  +  DMAG) 


iO  00 


c 


writeC12,996)  2,  2,  stf 
AUPPER(J)  =  AUPPER(J)  +  STF 
ALOWERCJ)  =  ALOWER(J)  +  STF 


C 

STF  =  RDAMP*(DDINC(1,N)*DDINC(2,N)) 
c  writeC12,996)  1,  2,  stf 

AUPPER(J-l)  =  AUPPER(J-l)  +  STF 
ALOWERCJ-1)  =  ALOWER(J-l)  +  STF 
220  CONTINUE 
300  CONTINUE 
C 

FORMATC'STFC ,I1,',',I1,')  =  *,E16.8) 
FORMATC'I IdeltaCdelta  D)ll  =  ’,E16.8) 
FORMATC'DOF  ',11,'  RHS  =  ',E16.8) 
FORMATC PARTICLE  NUMBER  ’,13) 

C 

RETURN 

END 


r>  o  o  o  o 


SUBROUTINE  ADDSTF(A , B , C , S , P , IDIAG , LM) 


ADDS  ELEMENT  STIFFNESS  AND  FORCE  TO  THE  GLOBAL  ARRAYS. 


C 

C 


20 

10 

C 


IMPLICIT  REAL*8(A-H,0-Z)  ^ 

REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 
DIMENSION  A(l) , B(l) , C(l) , IDIAG(l) 

DIMENSION  SC6,1),P(1),LM(1) 

DO  10  J=l,  6 

K=LMCJ)  • 

IFCK.EQ.0)  GOTO  10 

B(K)=B(K)+PCJ) 

L=IDIAG(K)-K 
DO  20  1=1,  6 
M=LMCI) 

IFCM.GT.K  .OR.  M.EQ.0)  GOTO  20 

M=L4M  • 

A(M)=ACM)+S(I,J) 

CCM)=C(M)+Sa,I) 

CONTINUE 

CONTINUE 


RETURN 

END 


c 


SUBROUTINE  NS0LVE(A , C , B,IDIAG , NEQ, FACT, BACK) 


....  PROGRAM  TO  PERFORM  (A/C)=L*D*U  FACTORIZATION  AND/OR 
BACKSUBSTITUTION  OF  AN  UNSYMMETRIC  SYSTEM  OF  EQUAHONS 
A(NA)  =  UPPER  TRIANGULAR  COEFFICIENT  MATRIX 
STORED  IN  COLUMN  FORM 

C(NA)  =  LOWER  TRIANGULAR  COEFFICIENT  MATRIX 

STORED  IN  COLUMN  FORM 

B(NEQ)  =  RIGHT  SIDE  VECTOR(AFTER  BACKSUBSTITION , 

IT  CONTAINS  THE  SOLUTION  VECTOR.) 

IDIAGCNEQ)  =  ADDRESSES  OF  DIAGONAL  TERM  IN  A(NA) 

NEQ  =  NUMBER  OF  EQUATIONS 

FACT  =  .TRUE.  ,  FACTOR  A  /  C 

.FALSE.  ,  DO  NOT  FACTOR  A  /  C 

BACK  =  .TRUE.  ,  FORWARD  REDUCE  B(NE0  AND  BACKSUBSTITUTE 
.FALSE.  ,  DO  NOT  FORWARD  REDUCE  B(NEQ)  OR 
BACKSUBTITUTE 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  DOT 

_  REMOVE  ABOVE  CARD  FOR  SINGLE-PRECISION  OPERATION 

LOGICAL  FACT, BACK 
DIMENSION  A(l) , B(l) , CCl) , IDIAG(l) 

_  FACTOR  A  TO  UT*D*U, REDUCE  B  TO  Y 

JR=0 

DO  10  J=1,NEQ 
JD=IDIAG(J) 

JH=JD-JR 

IF(JH.LE.l)  GOTO  10 
NS=J+1-JH 
NE=J-1 

I F(. NOT. FACT)  GOTO  20 
K=JR+1 
ID=0 

. . . .  REDUCE  ALL  EQUATIONS  EXCEPT  DIAGONAL 

DO  30  I=NS,NE 
IR=ID 

ID=IDIAG(I) 

NT=MIN0(ID-IR-1,I-NS) 

•  IF(NT.EQ.0)  GOTO  40 

A(K)=A(K)-  DOTCACK-NT),C(ID-NT),NT) 

C(K)=C(K)-  DOT(C(K-NT) , A(ID-NT) , NT) 

40  IF(A(ID).NE.0.)  C(K)=C(K)/A(ID) 

30  K=K+1 
C 


o  m  o  o  o  o 


C....  REDUCE  DIAGONAL  TERM  # 

C 

ACJD)=A(JD)-  DOT(A(JR+l),CCJR+l),JH-l) 

....  FORWARD  REDUCE  THE  R.H.S. 

20  IF(BACK)  B(J)=BCJ)-DOTCC(JR+l),B(NS),JH-l)  ^ 

10  JR=JD  • 

I F(. NOT. BACK)  RETURN 

....  BACKSUBSTITUTION 

J=NEQ 

JD=IDIAG(J)  • 

IF(A(JD).NE.0.)  B(J)=B(J)/A(JD) 

IFCNEQ.EQ.l)  RETURN 
50  D=B(J) 

J=J-1 

JR=IDIAG(J) 

IFQD-JR.LE.l)  GOTO  60 

M=J-JD+JR+2  • 

K=JR-M+1 
DO  70  I=M,J 
70  B(I)=BCI)-A(I+K)*D 
60  JD=JR 

IF(ACJD).NE.0.)  B(J)=B(J)/A(JD) 

IF(J.GT.l)  GOTO  50  • 

C 

RETURN 

END 


# 


o  o  nt  n  nt 


SUBROUTINE  MULT(A,B,C,IA,JA,JB) 

....  CALCULATE  THE  PRODUCT  C=A*B 
IMPLICIT  REAL*8CA-H,0-Z) 

....  REMOVE  ABOVE  CARD  FOR  SINGLE  PRECISION  OPERATION 

DIMENSION  A(IA,JA),  B(JA,JB),  C(IA,JB) 

C 

DO  5  1=1,  lA 
DO  10  J=l,  JB 
C(I,J)=0.0D0 
DO  15  K=l,  JA 

C(I,J)=Ca.J)+A(I,K)*B(K,J) 

15  CONTINUE 

10  CONTINUE 
5  CONTINUE 
RETURN 
END 


SUBROUTINE  MOVE(A,B,N) 

. . .  PROGRAM  TO  MOVE  N  ELEMENTS  OF  ARRAY  B  INTO  ARRAY  A 
IMPLICIT  REAL*8(A-H,0-Z) 

. . .  REMOVE  ABOVE  CARD  FOR  SINGLE-PRECISION  OPERATION 

DIMENSION  A(1),B(1) 

DO  10  1=1, N 
10  A(I)=B(I) 

RETURN 

END 


C  =========================:================= 

DOUBLE  PRECISION  FUNCHON  DOT(A,B,N) 

C 

C....  PROGRAM  TO  PERFORM  THE  DOT  PRODUCT  OF  TWO  VECTORS 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

C 

C _  REMOVE  ABOVE  CARD  FOR  SINGLE-PRECISION  OPERATION 

C 

DIMENSION  A(1),B(1) 

C 

DOT=0.D0 
DO  10  1=1, N 
10  DOT=DOT+A(I)*BCI) 

C 

RETURN 

END 


r>  o  r>  o  o 


c  - - - - -  0 

SUBROUTINE  CLEAR(A,M) 

....  PROGRAM  TO  CLEAR  A  FLOATING  POINT  ARRAY 


IMPLICIT  REAL*8(A-H,0-Z) 

. ..  REMOVE  ABOVE  CARD  FOR  SINGLE-PRECISION  OPERATION 

DIMENSION  ACl) 

DO  10  1=1, M 
10  ACI)=0.D0 

RETURN 

END 


SUBROUTINE  ICLEAR(IA,M) 

. . . .  PROGRAM  TO  CLEAR  AN  INTERGER  ARRAY 

DIMENSION  IA(1) 

DO  10  1=1, M 
10  IA(I)=0 

RETURN 

END 


